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FOREWORD 


This  final  report  describes  technical  work  accomplished  during  the 
Labyrinth  Seal  Analysis  program  conducted  under  Contract  F33615-80-C-2014. 
The  work  described  was  performed  during  the  period  June  1980  to  •9©“ 
W  85  .  This  contract  with  Allison  Gas  Turbine  Division  of  General 
Motors  Corporation  was  sponsored  by  the  Air  Force  Wright  Aeronautical 
Laboratories  Aeropropulsion  Laboratory,  United  States  Air  Force, 

Wright  Patterson  AFB,  Ohio,  with  Mr.  Charles  W.  Elrod  (AWAFL/POTX) 
as  Project  Engineer.  Technical  coordination  was  provided  by 
1st  Lt.  Keith  C.  Topham. 

The  technical  effort  reported  in  this  volume  was  performed  by 
personnel  of  Scientific  Research  Associates,  Inc.,  Glastonbury, 

Connecticut.  The  empirical  data  used  to  evaluate  the  results  of  the 
Analysis  Model  development  were  provided  by  Allison  Gas  Turbine  Division* 
This  report  was  submitted  in  four  volumes  in  May  1985.  Volume  I 
summarizes  the  development  of  the  labyrinth  seal  Analysis  Model.  Volume  II 
presents  the  user's  manual  for  the  Analysis  Model  computer  code. 

Volume  III  contains  the  experimental  results  and  summarizes  the  Design 
Model  based  on  these  empirical  data.  Volume  IV  presents  the  user's  manual 
for  the  Design  Model  computer  code* 

Publication  of  thla  report  does  not  constitute  Air  Force  approval  of 
the  findings  or  conclusions  presented.  It  la  published  only  for  the 
exchange  and  stimulation  of  ideas. 
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1.0  INTRODUCTION 


The  present  trend  of  aircraft  gas  turbine  design  has  been 
characterized  by  significant  increases  in  cycle  pressure  ratio  and  turbine 
inlet  temperatures  required  to  provide  higher  thermal  and  propulsive 
efficiencies.  Also,  increased  interest  in  engine  performance  and  fuel 
economy  has  created  additional  emphasis  for  improving  the  efficiency  of  gas 
turbine  engines.  These  trends  accentuate  the  need  for  improvements  In 
sealing  technology  and  the  development  of  advanced  design  and  analysis 
capabilities  to  reduce  gas  path  seal  leakage,  maintain  costly  vent  leakage 
to  a  minimum,  provide  better  control  over  sophisticated  cooling  circuits, 
and  prevent 'high  levels  of  seal  leakage  into  critical  aerodynamic  locations 
in  the  turbine  gas  path  which  can  result  in  considerable  penalty  from 
thermal  and  momentum  losses. 

Current  and  advanced  gas  turbine  engine  requirements  that  Impact 
labyrinth  seal  design  and  performance  include  a  broad  engine  power 
operating  range  which  usually  results  In  a  wide  range  of  seal  clearances* 
Normally,  seals  are  designed  to  run  as  tight  a  clearance  as  possible  at  the 
maximum  mission  time  condition.  In  setting  the  design  clearance, 
consideration  la  given  to  transient  differential  growth,  maneuver 
deflections,  mechanical  and  thermal  growths,  eccentricity,  and 
manufacturing  tolerances*  However,  with  variable  geometry  engines  and 
multiple  rolo  applications  the  engine  seals  will  not  always  operate  at  the 
design  clearance  nor  provide  minimum  leakage  across  the  operating 
spectrum.  Improved  seal  design  and  analysis  capabilities  must  be  developed 
to  address  this  problem. 

Cas  turbines  require  a  variety  of  labyrinth  seal  designs.  The  seal 
configuration  selected  for  a  given  application  is  based  on  the  purpose  of 
the  seal  and  satisfying  design  criteria  that  includes  the  following 
considerations:  axtal  envelope  available,  axial  travel,  clearance  range, 
potential  wear,  system  sensitivity  to  seal  clearances,  cooling  flow 
requirements,  sensitivity  to  damage  In  handling,  assembly  requirements,  and 
pressure  ratio. 
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Labyrinth  seals  are  used  throughout  a  gas  turbine  engine.  Including: 
compressor  and  turbine  airfoil  end  seals,  bearing  compartment  seals,  and 
flow  system  seals  to  minimize  or  control  flow.  The  purposes  of  these  seals 
are  not  always  the  same.  Labyrinth  seals  used  In  the  flow  path  are 
Intended  to  minimize  end  leakage.  Bearing  compartment  seals  are  Intended 
to  keep  the  oil  in  the  bearing  compartment  and  to  minimize  the  amount  of 
air  leakage  and  heat  addition  to  the  oil.  Thrust  balance  labyrinth  seals 
are  located  radially  to  provide  a  desired  off-setting  axial  load  component 
to  reduce  bearing  loads  to  the  design  level.  Other  flow  system  network 
seals  have  several  functions  Including:  controlling  leakage  flows  either 
to  a  minimum' or  to  a  level  to  satisfy  disc  pumping  and  thus  prevent  hot  gas 
recirculation  In  a  cavity,  controlling  cavity  pressures  to  reduce  axial 
bearing  loads,  or  preventing  excessive  leakage. 

The  variety  of  locations,  functions,  and  operating  conditions  Imposed 
on  labyrinth  seals  In  a  gas  turbine  engine  rt  quires  a  design  and  analysis 
capability  that  takes  advantage  of  the  numerous  seal  geometries  available 
and  accurately  predicts  the  seal  performance*  Labyrinth  seal  geometries 
Include  straight-through  seals,  step  seals ,  and  a  variety  of  advanced 
complex  geometry  designs.  The  seal  knives  may  be  vertical  or  slanted,  the 
knives  may  be  placed  on  a  rotating  or  stationary  surface.  Seal  lands  may 
be  smooth  and  solid,  honeycomb,  roughened  surface  solid,  striated,  r 
abradable  (porous  or  non-porous).  Other  seal  geosuetry  variables  Include 
knife  edge  thickness  and  sharpness,  clearance,  knife  pitch,  cavity  depth 
and  shape,  number  of  knives,  seep  height,  kntfe  location  on  the  step,  and 
knife  angle.  Aerodynamic  parameters  that  muat  be  considered  in  seal  design 
Include  rotational  speed,  pressure  ratio,  temperature,  and  Reynolds  number* 
U  a  labyrinth  seal  design  Is  to  be  successful  for  the  application 
Intended,  an  accurate  seal  design  and  analysis  model  Is  necessary.  The 
design  and  analysis  capabilities  available  today  re./  heavily  on  empirical 
relationships  which  severely  limit  the  application  range.  Available 
analytical  formulations  were  originated  many  years  ago  and  do  not  take 
advantage  of  modern  flow  field  calcuatlon  techniques  such  as  the  advances 
offered  by  solution  algorithms  for  the  Navtcr-Stokcs  equations.  Available 
seal  design  and  analysis  models  are  severely  restricted  relative  to 
analyzing  new  and  advanced  seal  designs.  Alan,  many  of  the  geometric  and 
aerodynamic  parameters  In  a  lab  seal  have  interfacing  effects  which  make  It 
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difficult  to  accurately  assess  Individual  parameter  effects  from  test 
data.  Therefore,  to  examine  the  numerous  Individual  and  combinations  of 
seal  geometric  and  aerodynamic  parameters  experimentally  would  be  time 
consuming  and  very  expensive.  In  addition,  empirically  derived  models  do 
not  provide  the  capability  to  assess  new  configurations  nor  do  they  provide 
the  design  engineer  with  guidance  on  how  to  improve  the  efficiency  of  the 
seal  beyond  what  information  has  been  determined  experimentally. 

Therefore,  a  critical  need  exists  for  labyrinth  seal  design  and 
analysis  calculation  models  that  provide  the  seal  design  specialist  with 
the  analytical  tools  to  calculate,  study,  understand  and  evaluate  the 
details  of  the  labyrinth  seal  Internal  flow  field  and  to  assess  subtle 
geometric  changes  relative  to  improving  the  seal  efficiency.  Final  tuning 
and  verification  of  the  resulting  configuration  would  still  be  accomplished 
on  a  seal  test  rig. 

Originally,  the  design  of  a  conventional  straight-through  labyrinth 
seal  was  usually  a  compromise  between  the  number  of  seal  knives  and  a  knife 
pitch  that  was  large  enough  to  reduce  the  kinetic  energy  carryover  to  a 
minimum.  However,  numerous  investigators  have  identified  a  significant 
number  of  additional  performance  Influence  parameters.  Today,  the 
qualified  seal  designer  recognises  there  are  a  large  number  of  geometric 
and  aerodynamic  parameters  that  Influence  the  performance  level  of  a 
labyrinth  seal.  These  parameters,  for  a  conventional  straight-through 
seal.  Include: 

o  Clearance  o  Pressure  ratio 

o  Pitch  o  Knife  angle 

o  Number  of  knives  o  Land  surface  (smooth,  honeycomb, 

o  Knife  tip  thickness  striated,  etc.) 

o  Rotational  speed  o  Knife  height 

o  Cavity  volume  o  Reynolds  number 

o  Knife  sharpness  o  Eccentricity 

To  the  above  list  several  additional  parameters  can  be  added  when  a 
step  seat  is  considered.  These  additional  parameters  Include: 

o  Step  height  o  Distance  fro*  seal  knife  to  step  face 

o  Step  configuration  o  Flow  direction 
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It  should  be  noted  that  the  seal  performance  Influence  parameters 
listed  above  do  not  operate  Independently.  The  Inherent  design  of  a 
labyrinth  seal  causes  the  Individual  geometric  and  aerodyamlc  Influence 
parameters  to  have  overlapping  or  Interfacing  effects.  For  example,  the 
effect  of  seal  clearance  on  leakage  Is  significantly  different  depending  on 
the  knife  pitch  and  number  of  knives.  Therefore,  a  very  large  matrix  of 
Identified  parameter  combinations  exists  which  determine  labyrinth  seal 
performance. 

There  are  numerous  types  and  applications  of  labyrinth  seals  in  a  gas 
turbine.  The  labyrinth  seal  may  be  a  straight-through,  slanted  straight,  a 
back-to-back  or  fir  tree  arrangement,  stepped,  or  slanted  stepped.  Flow 
may  be  "up”  or  "down”  the  step.  The  seal  lands  may  be  solid-smooth, 
roughened,  abradable  (porous  or  non-porous),  or  honeycomb. 

The  flow  fields  in  the  various  types  of  labyrinth  seals  used  in  a  gas 
turbine  have  similar  complexities  but  differ  significantly  between  a 
conventional  straight-through  and  stepped  seal  (see  Fig.  1).  The  stepped 
seal  has  a  mechanism,  the  vertical  step  face,  to  spoil  the  through  flow. 

The  straight  seal  has  a  core  of  through-flow  (referred  to  as  kinetic  energy 
carry-over)  that  results  in  higher  seal  leakage  rates  compared  to  a  stepped 
seal.  Therefore,  the  step  seal  provides  additional  parameters  for  the 
design  engineer  to  consider. 

The  number  of  tests  and  amount  of  hardware,  time,  and  cost  to  develop 
all  labyrinth  seal  performance  empirically  are  prohibitive.  In  addition, 
any  new  labyrinth  seal  design  concepts  would  have  to  be  tested  to  determine 
their  performance.  If  one  additional  geometric  or  aerodynamic  variable 
that  had  not  been  experimentally  evaluated  before  is  Introduced  or  varied, 
then  not  only  must  this  new  parameter  be  evaluated,  but  each  geometric  and 
aerodynamic  parameter  Interface  must  also  be  evaluated  to  determine 
combination  effects.  Therefore,  It  is  desirable  to  supplement  the 
empirical  approach  and  consider  analytical  techniques  to  assist  In  the 
deslgr  and  analysis  of  labyrinth  seals.  However,  the  labyrinth  seal  flow 
•field  must  be  classified  as  one  of  the  most  complex  and  challenging  for  a 
theoretical  analysis.  The  flow  field  Is  turbulent,  separated, 
compressible,  vIhcous,  has  one  wall  rotating,  and  experiences  streamwlse 
vorticlty.  If  unsteadiness  is  present,  then  any  analytical  analysis  that 
Is  not  time-dependent  may  have  difficulty  attaining  agreement  with  test 
results. 


4 


There  are  several  classical  analytical  methods  available  in  the 
literature  for  estimating  labyrinth  seal  leakage.  However,  erch  of  these 
methods  is  based  on  certain  simplifying  assumptions  which  limit  the  areas 
of  applicability.  Several  of  the  more  recent  theoretical  methods  will 
yield  reasonable  estimates  of  leakage,  but  the  area  of  applicability  of 
each  method  is  restricted  to  a  narrow  range  of  geometries  and  overall 
pressure.  None  of  the  available  methods  account  for  more  than  four  or  five 
variables.  No  general  solution  is  available  for  calculating  labyrinth  seal 
leakage  flows,  nor  has  an  analytical  method  been  developed  that  examines 
the  seal  interior  flow  field  to  provide  aerodynamic  details  and  guidance  to 
improve  the  efficiency  of  the  labyrinth  seal. 

The  origin  of  the  labyrinth  seal  can  be  traced  to  C.A.  Parsons 
(Ref.  1)  who  apparently  introduced  the  concept  of  a  steam  turbine  and 
reported  the  event  in  1892  (Ref.  2).  The  design  concept  was  to  provide  a 
tortuous  path  between  high  and  low  pressure  regions  by  using  a  series  of 
non-contacting  restrictions  and  chambers.  The  characteristic  functions  of 
the  restrictions  were  to  convert  the  pressure  head  into  kinetic  energy 
which  would  then  be  dissipated  as  completely  as  possible  in  the  intervening 
chambers.  The  effecti1 2  ne«s  of  this  concept  is  shown  by  its  continued  use 
in  the  current  most  mouern  gas  turbine  designs. 

Although  the  labyrinth  seal  is  relatively  simple  in  design,  the 
numerous  geometric  and  aerodynamic  parameters  associated  with  determining 
the  overall  performance  are  numerous  and  complex  as  noted  earlier. 

However,  theoretical  formulations  to  describe  the  flow  field  through  the 
labyrinth  seal  have  been  attempted  apparently  starting  with  E.  Becker  in 
1907  (Ref.  3).  Closely  following  Becker  was  a  paper  by  H.M.  Hartin 
(kef.  4).  It  is  interesting  to  note  that  these  two  papers  establish  the 
classical  labyrinth  seal  theories  which  can  be  organized  into  two  families; 

(1)  Treat  the  seal  knives  as  a  series  of  individual  throttles 
(Hartin  approach,  Ref.  4). 

(2)  Treat  the  seal  as  a  friction  device  (Becker  approach.  Ref.  3). 
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Generally  speaking,  existing  labyrinth  seal  flow  calculation  theories 
make  assumptions  that  place  the  theories  in  one  of  two  classifications 

(1)  Small  AP  is  assumed  across  each  restriction. 

(2)  Last  restriction  choked. 

Stodola's  (Ref.  5)  and  Martin's  (Ref.  4)  formulas  assume  small  AP  across 
each  restriction,  i.e.,  all  kinetic  energy  is  recovered  as  Internal  heat  in 
the  seal  cavities.  These  formulas  also  assume  the  discharge  coefficient  is 
1.0  for  each  restriction.  These  two  often-quoted  formulas,  as  well  as 
others  that  are  similar,  assume  that  the  gas  experiences  an  ideal 
isothermal  expansion  across  each  seal  knife  followed  by  dispersion  of  the 
kinetic  energy  and  reheating  before  entering  the  next  seal  knife.  The  many 
theoretical  studies  in  this  classification  which  have  appeared  in  the 
literature  are  of  little  value  in  the  prediction  of  leakage  rates  through 
the  commonly  used  straight  labyrinth  seal  because  the  assumptions  on  which 
they  are  based  approach  the  conditions  which  are  approximated  only  in  step 
seals.  The  authors  of  these  theoretical  studies  seek  to  modify  the  theory 
with  empirical  correction  factors  to  make  the  calculation  fit  test  data. 

Theoretical  analyses  to  date  have  neglected  or  approximated  the 
carryover  (sometimes  referred  to  as  velocity  of  approach)  of  kinetic  energy 
from  knife  to  knife.  This  carryover  factor  varies  substantially  with 
geometry  and  pressure  ratio.  Also,  the  discharge  coefficient  must  be 
evaluated  for  each  particular  knife  (or  stage)  geometry  and  will  vary 
depending  on  several  factors  including: 

o  Knife  tip  thickness  o  Reynolds  number 

o  Clearance  o  Number  of  Knives 

In  addition  to  the  foregoing  parameters,  the  carryover  factor  and 
discharge  coefficient  may  be  effected  by  surface  roughness,  land  porosity, 
honeycomb  and  striated  lands,  and  knife  rotation. 


Trutnovaky  (Ref.  6).  like  Becker  (Ref.  3).  treated  the  seal  leakage 
as  flow  in  a  rough  pipe.  However,  his  solution  of  the  basic  equations 
describing  the  flow  is  complicated  and  in  general  difficult  to  use. 
Zabri8kie  and  Sternlicht  (Ref.  7),  offer  an  approach  that  is  an  extension 
and  simplification  of  the  method  used  by  Trutncvsky.  However,  numerous 
investigators  have  challenged  this  general  approach  as  not  being  correct 
relative  to  the  physical  considerations  of  the  problem. 

Most  theoretical  approaches  make  a  simplifying  assumption  regarding 
tho  intercavity  vortex  and  eddies  and  give  no  help  on  how  to  improve  these 
cavitlen  to  reduce  leakage.  It  is  known  that  the  shape  and  size  of  the 
cavity  between  seal  knives  affects  the  strength  of  the  vortex  and  eddies 
which  convert  the  kinetic  energy  issuing  from  each  knife  into  internal 
energy. 

It  is  spp^ren-  from  the  foregoing  discussion  that  existing 
theoretical  methods,  which  necessarily  employ  empirically  derived  modifiers 
for  calculating  labyrinth  seal  flows,  have  significant  shortcomings.  A 
generalized  theoretical  approach  with  proven  accuracy  and  reliability  is 
not  available. 

This  situation  leaves  the  design  engineer  with  the  dubious  task  of 
selecting  a  method  for  hie  needs  from  numerous  methods  with  varying  degrees 
of  accuracy  and  range  of  application.  Furthermcre,  the  seal  design 
specialist  does  not  Have  the  analytical  tools  to  examine  the  details  of  the 
labyrinth  seal  interior  flow  field  ana  determine  geometry  changes  that 
would  increase  seal  cavity  turbulence  resulting  in  improved  sealing 
efficiency. 

The  analysis  tools  required  by  the  mechanical  design  engineer  and  the 
sea)  research  specialist  are  distinctly  different.  The  design  engineer 
requires  a  simplified  "design"  calculation  model  that  will  determine  the 
overall  performance  of  a  labyrinth  seal  when  selected  geometry  la 
specified.  The  design  engineer  also  has  a  need  for  a  calculation  aodel 
that  will  provide  dimensional  criteria  for  an  optimum  seal  configuration 
for  a  given  application  when  an  absolute  minima,  of  information  is 
supplied,  i.e.,  clearance,  axial  envelope,  otatlonal  speed,  and  pressure 
ratio.  The  design  model  input  format  should  be  simple  and  the  computer  run 
time  should  be  minimum  for  the  model  to  be  practical  as  a  production 
program.  Incorporating  available  test  data  and  generating  seal  nar romance 
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data  not  currently  available  to  update  and  expand  existing  theoretical 
models  that  rely  heavily  on  empirical  correlations  offers  the  most 
effective  approach  for  the  development  of  an  advanced  "design"  model. 

The  seal  research  specialist  requires  an  “analysis"  calculation  model 
that  provides  the  aerodynamic  details  of  the  seal  interior  flow  field  in 
order  to  determine  and  evaluate  the  effects  of  the  numerous  geometric  and 
aerodynamic  parameters  incorporated  in  the  design  of  conventional 
straight-through  and  stepped  seals.  This  capability  will  enable  the  seal 
specialist  to  identify  and  analytically  evaluate  design  improvements  to 
obtain  higher  efficiency  labyrinth  seals  as  well  as  to  improve  the  accuracy 
of  conventional  seal  design  leakage  calculations. 

Recent  developments  in  the  phenomenological  models  of  internal 
turbulent  flows  shows  that  the  methodology  has  progressed  to  the  point  that 
the  complex  turbulent  flow  field  within  the  labyrinth  seal  interior  may  be 
calculated  via  direct  analysis  using  Navier-Stokes  computer  codes  presently 
available,  modified  for  the  geometry  of  a  labyrinth  seal.  The  successful 
application  of  a  compressible  time  dependent  Navier-Stokes  solution  would 
provide  a  major  breakthrough  in  seal  analysis  technology. 

The  value  of  a  Navier-Stokes  solution  method  for  the  labyrinth  seal 
leakage  calculation  is  that  it  can  potentially  analyse  most  of  the 
geometric  and  aerodynamic  effects  individually  and  in  matrix  combinations. 
There  may  be  some  effects  that  cannot  be  completely  modeled.  In  these 
cases,  test  data  may  be  used  to  support  and  expand  the  Navier-Stokes 
solution  method. 

Labyrinth  seal  design  improvements  have  been  limited  because  the 
tools  to  analyze  the  effects  of  geometric  changes  or  unique  configurations 
do  not  exist  except  In  a  very  fundamental  or  simplified  analysis  form.  The 
availability  of  a  Navier-Stokes  solution  would  provide  a  capability  to 
study  and  analyze  many  complex  geometric  shapes  and  configurations  that  can 
only  be  evaluated  presently  through  expensive  and  time  consuming  tests. 
Although  the  Navier-Stokes  solution  may  be  limited  in  calculating  the  exact 
level  of  performance  for  exotic  seal  configurations,  it  will  guide  the 
engineer,  through  performance  trends,  to  a  more  efficient  design.  The 
final  design  should  be  tested  to  verify  performance  characteristics  of  the 
seal . 
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In  regard  to  Navier-Stokes  solutions,  several  items  need  to  be 
considered*  It  is  clear  that  a  viable  analysis  which  simulates  the  seal 
flow  field  must  allow  for  flow  ranging  from  the  low  subsonic  regime  to  the 
transonic  regime,  must  Include  possible  shock  waves,  must  allow  for 
dominating  viscous  effects  and  must  allow  for  very  large  regions  of 
recirculation.  These  considerations  clearly  dictate  a  Navier-Stokes 
approach  to  the  problem* 

In  most  instances  the  Navier-Stokes  equations  are  so  intractable  that 
only  numerical  solutions  can  be  obtained.  Numerical  techniques  for  solving 
the  Navier-Stokes  equations  are  discussed  by  Roache  (Ref.  8)  and  more 
particularly  for  the  compressible  Navier-Stokes  equations  by  Peyret  and 
Viviand  (Ref.  9).  Peyret  and  Viviand  singled  out  three  techniques,  the 
explicit  scheme  of  MacCormack  (Ref.  10),  the  scheme  due  to  Wldhoff  and 
Victoria,  (Ref.  11)  also  explicit,  and  the  implicit  scheme  of  Briley  and 
McDonald  (Ref.  12).  The  technique  of  MacCormack  (Ref.  10)  has  been  very 
effectively  applied  for  instance  by  Shang,  Hankey  and  Law  (Ref.  13)  in  a 
shock  wave-boundary  layer  Interaction  problem.  However,  the  need  to 
compute  large  regions  of  relatively  low  speed  recirculating  fluid 
interspersed  with  local  high  speed  throats  in  the  labyrinth  seal  problem 
could  make  the  stability  bounds  of  the  MacCormack  scheme  very  restrictive. 

As  a  result  in  the  labyrinth  seal  problem,  unacceptably  long  computer 
run  times  could  result  from  the  required  locally  refined  spatial  meshes 
with  the  stability  restricted  scheme.  The  technique  of  Briley  and 
McDonald  (Ref.  12)  on  the  other  hand,  is  not  restricted  by  the  stability 
bounds  of  the  MacCormack  scheme  and  the  Wldhoff  and  Victoria  scheme  and, 
thus,  is  better  suited  for  the  labyrinth  seal  problem  than  either  of  the 
other  two  candidate  algorithms.  For  these  reasons  the  Briley-McDonald 
(Ref.  12)  technique  was  used  In  the  present  effort  to  predict  the  detailed 
flow  field  In  the  labyrinth  seal  investigation. 
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2.0  ANALYSIS 


2.1  Governing  Equations 

The  governing  equations  utilized  in  this  study  are  the  ensemble 
averaged  time  dependent  Navier-Stokes  equations.  These  equations  are  the 
mathematical  statement  of  the  physical  conservation  laws  of  mass,  momentum 
and  energy  for  one  phase  fluid  dynamic  systems.  Using  vector  notation, 
these  transport  equations  can  be  respectively  written  as 

dp  _ 

—  *  -V  (pV)  (1) 

a  „  ™ 

-  (f  V)  =  -  V(/>W)  +  V- t  -  VP  (2) 


and 


~  V-  V)]  »  -7-  [/>(U  +  T>  H)v]  -  V  q  +  V*(r  *  V)  -  7  (PV) 


(3) 


♦  * 
were  p  is  the  fluid  density,  V  Is  the  Eulorlan  velocity  vector  and  U  la 

the  specific  Internal  energy.  The  thermodynamic  pressure,  P,  ia  related  to 

the  temperature,  T,  and  density,  p,  by  the  perfect  gas  equation  of  state 


P  * 


R 

*>irT 


(«> 


where  R  Is  the  universal  gas  constant  and  M  Is  the  mixture  molecular 
weight.  The  stress  tensor,  t.  Is  related  to  the  velocity  vector,  by  the 
relationship 


2,  _  — 


r  »/i{7V4-Vvt  - 


(5) 
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where  w  is  the  viscosity  while  the  heat  transfer,  q,  is  related  to  the 
temperature,  T,  by  the  Fourier  relationship 

q’  =  -«VT  (6) 

where  <  is  the  thermal  conductivity.  Eq.  (3),  the  energy  equation, 

represents  the  balance  between  the  time  rate  of  charge  of  the  Internal  plus 
*  ♦  ♦ 

kinetic  energy  (U+l/2V«V),  the  convection  of  that  energy,  the  heat  transfer 
and  the  stress  and  pressure  work.  Defining  the  static  enthalpy  by 

a  p 

h  a  U  +  —  (7) 

P 

the  energy  equation  can  be  rewritten  as 


+  ^  V  V)]  -  •—  =  -V-[/>(h  +  ~V-  V)]-  V-q*  +  V-(rV) 


(8) 


Using  the  stagnation  enthalpy  defined  by  the  relationship 


h0  «  h  +~V-V 


(9) 


yields  the  more  compact  form,  viz.. 


dt 


(ph) 


dP 

IT 


*  -V'(p  h0V )  -  +  V*  (  t’ V) 


(10) 


By  dotting  the  velocity  vector,  0,  with  the  momentum  equation,  Eq.  (2)  and 
using  the  vector  Identities 


—  d  ~  d  ,  V  V  ,  V  V  dP 

v  3T 


(ID 


and 


-  — .  r  V  •  V  —-i  V  V  - 

VV(pW)  =  V •[(/>“ )VJ  ♦  — 2“  V  {pV)  (12) 
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and  applying  the  continuity  equation,  Eq.  (1),  one  can  obtain  the  so-called 
mechanical  energy  equation 


d  VV  r  V*V  —i  —  - 

—  (/) — —  [(/>  — ~)VJ- V-VP  +  V  (V  r) 


(13) 


Subtracting  the  mechanical  energy  equation,  Eq.  (13),  from  the  energy 
equation,  Eq.  (8),  and  applying  the  vector  identity 

V*(r  V) s  V  (V  t)  +  t  :  VV 
and  using  the  definition  of  energy  dissipation,  ♦  ,  as 

$  s  r  t  V V 

one  obtains  an  alternate  form  of  the  energy  equation  expressed  in  terns  of 
the  time  rate  of  charge  of  the  static  enthalpy 


(14) 


(15) 


a  dP  _  _ 

dT(^h)"  V*q  +  V'7P  +  <J> 


(16) 


The  static  enthalpy,  h,  can  be  related  to  the  temperature,  T,  by  the 
relat ionship 


T 

h  *  /  Cp(T>  dT 

Tr 


(17) 


where  Cp  is  the  specific  heat  at  constant  pressure.  If  a  calorically 
ideal  gas  is  assumed,  (cp  is  a  constant). 


h  »  CpT 


(18) 


Eq.  (16)  can  be  further  simplified  and  written  in  terms  of  temperature. 
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The  system  of  coupled  nonlinear  partial  differential  equations 
represented  by  Eq.  (1),  (2),  and  (16)  are  the  basis  of  the  governing 
equations  used  in  this  study.  The  equations  are  valid  for  both  laminar  and 
turbulent  flow.  However,  for  turbulent  flow,  all  variables  are  ensemble 
averaged  (Refo  14)  and  the  viscosity  and  thermal  conductivity,  k,  must  be 
considered  as  effective  values.  Thus,  the  viscosity  must  be  considered  as 
the  sum  of  the  laminar  and  turbulent  viscosities  (the  turbulent  viscosity 
comes  from  the  Boussinesq  approximation  to  the  Reynolds  stress  terms). 


M  *  /*,  + 


(19) 


The  viscosity  is  related  to  the  thermal  conductivity  and  the  specific  heat 
by  the  concept  of  a  Prandtl  number  (which  is  presumed  to  be  known) 


ki  = 


CP^i 


Pr. 


kT  * 


CP^T 


Pr, 


and 


(20) 


(21) 


k  *  k|  +  kT  (22) 

In  this  study,  two  methods  were  used  to  represent  the  turbulent  viscosity, 
yt.  The  first  uses  the  algebraic  mixing  length  of  Prandtl  whore 

fir  83  pt m  o/OiD  (23) 


where  D:t)  is  the  second  variant  of  the  mean  flow  rate  of  deformation 
tensor,  i.o.. 


0:0 


$ 


♦  *§■  (7- v  )2 


(24) 


and  tm  is  the  mixing  length  which  will  be  discussed  at  a  later  time. 
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The  second  method  utilized  in  this  study  is  to  assume  the 
Prandtl-Kolmogorov  relationship  (Ref.  15)  for  turbulent  viscosity,  the 
so-called  k-e  model 

H-j  -C^—-  (25) 


where  is  a  'constant'  (to  be  discussed  at  a  later  time)  k  is  the 
turbulence  kinetic  energy  defined  by 


(26) 


and  e  is  the  dissipation  of  turbulence  kinetic  energy.  Partial 
differential  equations  govern  the  distribution  of  k  and  e  in  the  flow.  For 
a  discussion  of  these  equations,  the  reader  is  referred  to  Refs.  16 
and  17.  The  resultant  equations  are 


dpK 

~aT~ 


«  -V-(^VK)+V-(~V K)+  (/iT 


0;  0 


-/>«) 


(27) 


Opt  LLr  t 

-  *  K)+V*  ( — _7«  )  ♦  — (C,/1T0:0  -  C tpt)  (28) 

it  cr(  K 

Thus,  when  a  k-e  turbulence  model  Is  used,  two  additional  coupled  nonlinear 
partial  differential  equations  are  added  to  the  original  governing  partial 
differential  equations. 

In  summary,  the  governing  system  of  equations  used  In  this  study 
include  the  partial  differential  equations,  Kqs.  (1),  (2),  and  (16)  plus 
(27)  and  (28)  If  a  k-e  turbulence  model  is  used.  The  constitutive 
relationships  are  represented  by  Kqs.  (4)-(7),  (15),  (17)  or  Kqs. (18), 
(l9)-(22),  and  (23)  If  a  mixing  length  model  Is  used,  or  Kq.  (25)  If  a  k-e 
turbulence  model  is  used. 
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2.2  Coordinate  Systems 


Application  of  the  governing  system  of  partial  differential  equations 
to  a  given  problem  is  not  in  general  straight  forward  and  several  decisions 
must  be  made  before  these  equations  can  be  put  in  a  form  suitable  for 
solution  by  numerical  techniques.  First,  a  coordinate  system  must  be 
chosen.  For  the  cases  of  interest  in  this  study,  a  body  or  boundary 
conforming  coordinate  system  is  normally  used.  Except  for  the  simplest 
cases,  a  noncartesian  coordinate  system  must  be  employed.  In  many  cases  a 
nonorthogonal  coordinate  system  must  be  used,  e.g.,  the  case  of  a  seal  with 
a  tapered  knife  blade.  Use  of  a  noncartesian  coordinate  system  requires 
the  choice  of  both  the  components  of  velocity  vector,  V,  and  the  choice  of 
the  directions  in  which  the  vector  momentum  equation  is  to  be  expressed* 

For  instance,  the  velocity  components  and  the  directions  in  which  the 
momentum  equations  are  written  can  be  aligned  with  the  coordinate 
directions  (Ref.  18),  or  if  one  desires,  the  velocity  components  and 
momentum  equation  directions  could  be  aligned  with  the  cartesian  or 
cylindrical  polar  coordinate  directions  (Refs.  19-21).  Combinations  of  the 
above  are  also  possible,  although  to  the  author's  knowledge,  have  not  been 
used  to  date. 

In  this  study  two  basic  types  of  geometric  configurations  are 
considered:  (1)  planar  configurations  and  (2)  axiaysmetrlc 
configurations.  The  planar  configurations  arc  utilized  primarily  to 
simulate  experimental  setups  where  the  data  Is  taken  In  a  two-dimensional 
planar  environment.  The  axtsymmetrlc  configurations  are  used  primarily  to 
simulate  engine  component  performance  tn  a  rotating  or  nonrotating  (but 
still  axtsymmetrlc)  environment.  When  a  rotating  system  Is  analyzed,  the 
flow  Is  three-dimensional,  however,  there  is  symmetry  with  respect  to  the 
rotating  direction,  t.e.,  2/39-0.  Thus  the  governing  partial  differential 
equations  can  be  expressed  in  terms  of  two  coordinate  directions,  but  three 
velocity  components  are  required  to  define  the  velocity  vector  and  three 
momenta  equations  must  be  solved. 

The  restalnder  of  this  section  will  describe  the  derivation  of  the 
governing  system  of  partial  differential  equations  to  he  used  for  planar 
nonorthogonal  configurations.  The  derivation  for  axtsymmetrlc 
configurations  Is  similar  and  the  details  arc  discussed  In  Ref.  19. 
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Although  in  this  study  the  governing  equations  were  only  solved  in  two 
dimensions  and  axisymmetric  flow,  the  technique  is  applicable  to  two  and 
three  space  dimensions.  For  generality  here  the  three-dimensional 
derivation  is  presented. 

The  technique  utilized  in  this  study  is  similar  to  that  reported  in 
Refs.  19-21.  The  dependent  variables  are  chosen  as  the  three  cartesian 
velocity  components,  uj ,  u2  and  U3,  the  density,  p,  and  the  static 
enthalpy,  h,.  The  momentum  equation  is  also  expressed  in  terras  of  the 
three  cartesian  directions  xj,  x2  and  xj.  These  governing  partial 
differential  equations  in  cartesian  coodinates  can  be  expressed  in  the  form 
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where  in  tensor  notation,  the  nonaal  and  shear  stress  component#  and  the 
heat  transfer  are  respectively 
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and 


(39) 


q*i  dxj 

Transforming  Eq.  (30)  to  a  "caeral  coordinate  system  in  which  the 
general  coordinates,  y J ,  are  related  to  the  cartesian  coordinates 
xi,  X2  and  X3  by  the  relationship 

yJ  s  yj< vvx3)  v 

and  use  of  the  chain  rule  yields 

aw  ap,  ay' 

=  (l 

Defining  the  Jacobian  determinant  of  the  inverse  transformation,  J,  by 
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and  multiplying  Eq.  N41)  by  the  Jacobian  yields 


dW  dy  dF, 

J—  =  J  — - 7  (A3) 

dt  dx(  dyj 

This  form  is  sometimes  referred  to  as  the  weak  conservative  form.  By 
assuming  that  the  Jacobian  is  not  a  fir.-  tion  of  time,  t,  one  can  rearrange 
Eq.  (AA)  in  the  form 

d  .  d  ,  ay1  „  t  ..  d  ,  ,  dyJ  . 


—  (JW)  *  — 
dt  dy 


)  dxj  dyi  dx, 


18 


However,  It  can  be  shown  that  (Appendix  A) 


-^y  (J  ill.)  =  0  (45) 

ayJ  dXj 


The  derivation  of  Eq.  (44)  in  fact  is  general  in  that  if  the 
transformation  of  Eq.  (43)  had  a  moving  coordinate  system,  the  Jacobian 
could  be  taken  inside  the  time  derivative,  and  thus  the  transformed 
equation  can  be  written  in  the  form 


d 

*r~ ( JW  #  * 

dt 


(46) 
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This  form  of  the  transformed  governing  partial  differential  equations  is 
sometimes  referred  to  as  the  strong  conservative  form,  and  this  is  the  form 
of  the  governing  equations  solved  In  this  study.  The  primary  advantage  of 
using  the  cartesian  velocity  components  and  solving  the  scalar  momenta 
equations  associated  with  the  cartesian  directions  is  that  the  number  of 
terms  of  the  governing  differential  equations  Is  kept  to  a  minimum  for 
nonorthogonal  systems.  If  there  is  no  specific  reason  why  the  carte. tan 
velocity  components  and  directions  for  the  momenta  equations  are  unsuitable 
for  a  given  application,  this  Is  the  most  efficient  means  of  solving  the 
governing  equations.  For  this  study  and  many  other  applications  this 
technique  does  not  appear  to  have  any  disadvantages  vis  a  vis  other 
methods,  and  hence  that  procedure  was  chosen  here. 

In  three  dimensions  each  convective  and  pressure  term  from  the 
cartesian  equations  In  general  becomes  three  terms  In  the  transformed 
system.  Each  stress  and  heat  transfer  term  becomes  nine  terms  in  the 
transformed  system.  For  two-dimensional  flow,  the  corresponding  numbers 
are,  two  and  four*  To  define  the  stress  and  heat  transfer  terms  in  the 
transformed  coordinate  system  requires  applying  the  transformation  to 
Eqs.  (35-39).  This  yields  In  tensor  notation, 
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and 


(50) 


2.3  Initial  and  Boundary  Conditions 

Steady  solution  of  the  system  of  governing  partial  differential 
equations  represented  by  Eq.  (46)  is  obtained  by  time  marching  this  system 
of  equations  until  the  steady  state  is  reached.  Before  the  solution 
procedure  is  described  two  important  aspects  must  be  discussed:  (1)  the 
initial  conditions  and  (2)  the  boundary  conditions.  Any  procedure  which 
utilizes  either  a  time  marching  method  to  obtain  a  steady  state  (or 
transient)  solution  or  a  Newton-Raphson  Iteration  procedure  requires  some 
initial  guess  of  the  flow  variables  (in  this  case  all  the  dependent 
variables  and  other  necessary  variables  such  as  pressure,  temperature, 
viscosity,  etc.).  In  some  of  the  simpler  cases,  a  reasonable  approximation 
to  a  converged  solution  can  either  be  guessed  or  obtained  through  physical 
reasoning.  However,  since  the  range  of  geometries  considered  under  this 
effort  were  so  diverse,  It  was  felt  that  no  reasonable  general  Initial 
guess  procedure  could  be  developed.  The  approach  taken  here  was  to  assure 
that  the  flow  was  Initially  stagnant  (all  velocity  components  were  zero) 
by  assuming  that  the  pressure  and  temperature  were  a  constant  and  set  equal 
to  the  stagnation  conditions  of  the  source  flow.  The  downstream  or  back 
pressure  is  then  lowered  to  some  desired  level  over  a  period  of  time,  and 
the  flow  Is  drav.ii  through  the  seal  until  a  steady  state  is  achieved. 
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This  technique  has  the  advantage  of  being  easy  to  implement  in  any 
geometric  configuration.  In  addition,  the  lowering  of  the  back  pressure 
can  be  considered  to  be  similar  to  an  experimental  apparatus  where  the 
source  tank  is  pumped  up  and  the  back  pressure  gradually  lowered  by  opening 
a  downstream  valve. 

To  obtain  a  solution  of  the  governing  system  of  partial  differential 
equations  represented  by  Eq.  (46),  it  is  necessary  to  define  boundary 
conditions  on  each  bounding  surface  of  the  computational  domain.  For  the 
purposes  of  -this  investigation  boundary  conditions  can  be  classified  as 
occurring  on  three  different  types  of  bounding  surfaces:  (1)  walls  or 
solid  surfaces,  (2)  inlets  and  (3)  exits.  The  boundary  conditions  utilized 
on  each  different  type  of  surface  will  now  be  discussed  in  turn.  Analysis 
of  the  characteristics  of  the  boundary  layer  equations  shows  that  (in  three 
dimensions)  four  conditions  must  be  specified  on  walls.  For  this  study, 
the  no-slip  conditions  are  used  for  the  tangential  (to  the  wall)  velocity 
components,  i.e., 


»  0 


(51) 


and 


U 
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(52) 


where  the  subscripts  Tj  and  T2  refer  to  the  two  tangential  velocity 
components.  (For  two-dimens  1 on a  1  flow  only  one  tangential  velocity 
component  is  used).  For  the  normal  velocity  component  either  the  normal 
velocity  component  or  the  normal  mass  flux  is  specified,  i.e., 
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where  the  subscript  w  refers  to  the  specified  wall  value.  The  fourth 
condition  used,  the  thermal  condition,  either  specifies  the  wall  as  being 
adiabatic  or  specifies  che  wall  temperature.  These  conditions  can  be 
written  respectively  as 

rTw  •  V T  =  0  (55) 

or 

T  *  Tw  (56) 

where  in  this  case  nw  represents  the  unit  vector  normal  to  the  wall.  In 
addition,  a  fifth  condition  (in  two-dimensions  a  fourth)  not  required  by 
the  characteristic  analysis,  is  used  for  convenience  to  close  the  set  of 
equations.  The  need  for  this  fifth  condition  could  be  removed  by  the  use 
of  one-sided  difference  approximations  or  by  applying  one  of  the  governing 
equations  at  the  wall.  In  this  study,  the  second  method  was  used  and  the 
boundary  layer  approximation  to  the  normal  momentum  equation  was  applied  at 
the  wall.  This  can  be  expressed  as 

rw-VP«0  (57) 

Studies  have  Indicated  that  there  is  little  difference  between  using  this 
equation  or  using  the  full  normal  momentum  equation.  The  condition  of 
Eq.(57)  approximates  the  normal  momentum  equation  to  order  Re”1  for 
viscous  flow  at  a  no-sllp  surface.  The  symmetry  equations  are  meant  to  be 
applied  on  a  plane  or  axis  of  symmetry.  The  normal  component  of  velocity 
Is  set  equal  to  zero  on  the  axis  of  symmetry,  l.e., 

(TS  V*0  (58) 

♦ 

where  n8  is  the  unit  vector  normal  to  the  axis  or  plane  of  symmetry*  In 
addition  the  normal  derivatives  of  the  remaining  two  tangential  components 
of  velocity  arc  set  to  zero.  Two  other  conditions  must  be  set  on  the  axis 
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or  plane  of  symmetry.  Usually  the  symmetry  conditions  on  pressure  and 
temperature  are  used,  viz., 


n$  VP«0  (59) 

and 

n#  *  VT  «  0  (60) 

On  inlets  the  characteristic  analysis  shows  that  four  conditions  must 
be  set  (for  three-dimensional  flow).  The  procedure  used  in  this  study  is 
to  divide  the  flow  on  the  inlet  into  two  regions:  (1)  a  central  core  where 
the  flow  is  essentially  inviscid  and  (2)  the  attached  boundary  layers  where 
the  normal  pressure  gradient  is  zero.  This  technique,  which  can  be  called 
the  ’two-layer  model’,  sets  a  constant  stagnation  pressure  in  the  central 
core  region  of  the  inlet  boundary.  In  the  attached  boundary  layer(s),  the 
static  pressure  is  set  at  the  central  core  edge  value,  and  the  form  of  the 
streamwise  velocity  profile  is  set.  It  is  to  be  emphasized  that  the 
magnitude  of  the  streamwise  velocity  is  not  specifically  set.  Rather 
interaction  of  the  core  flow  and  the  rest  of  the  flow  with  the  inlet 
boundary  layers  sets  the  magnitude  of  the  boundary  layer  inlet  streamwise 
velocity  profile.  In  this  study  two  forms  of  boundary  layer  profiles  were 
used:  (1)  for  laminar  flow  the  von-Karman-Polhausen  (Ref.  22)  profile  was 
used,  and  (2)  for  turbulent  flow  the  Maise-McDonald  (Ref.  23)  profile  was 
used.  In  addition  to  the  ’two  layer  model'  boundary  condition,  two 
additional  velocity  conditions  must  be  specified.  In  this  study  the  flow 
angle  between  the  streamwise  velocity  component  and  the  transverse 
component  was  specified  (usually  as  zero)  and  if  the  apparatus  was  rotating 
a  swirl  velocity  profile  was  specified*  The  fourth  boundary  condition  was 
the  thermal  condition  that  the  stagnation  enthalpy  on  the  Inlet  boundary 
remains  constant.  In  addition,  a  fifth  condition,  not  required  by  the 
characteristic  analysis.  Is  used  to  numerically  close  the  set  of 
equations.  In  this  case  the  weak  condition  that  second  derivative  of 
pressure  normal  to  the  inlet  plane  equals  zero  was  used.  The  advantage  of 
this  condition  Is  that  It  allows  pressure  waves  to  exit  the  computational 
domain  rather  than  reflect  off  boundaries. 
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A  characteristic  analysis  of  the  governing  equations  on  an  exit  plane 
shows  that  for  subsonic  flow  only  one  condition  must  be  set.  In  this  study 
the  condition  is  met  by  setting  the  back  pressure  to  some  desired  value. 
Thus,  since  the  upstream  (inlet)  stagnation  pressure  is  set  in  the  inlet 
boundary  core  region,  the  mass  flux  is  determined  by  these  two  variables 
and  the  loss  mechanism  that  occurs  in  the  physical  domain.  As  before  to 
numerically  close  the  set  of  equations,  four  more  exit  conditions  mus'-  be 
set.  In  this  case  weak  conditions  are  set,  viz.,  the  second  derivative  of 
the  three  velocity  components  and  the  temperature  are  set  to  zero,  the 
so-called  parabolic  assumptions* 

If  the  two-equation  (k-e)  turbulence  model  is  used,  two  additional 
partial  differential  equations  must  be  solved,  and  hence  initial  and 
boundary  conditions  must  be  specified  for  these  equations.  The  procedure 
used  in  this  study  is  to  obtain  a  converged  solution  with  a  mixing  length 
model,  then  use  the  assumption  of  an  'equilibrium  turbulence  model*  as  an 
initial  guess  for  the  k  and  e  fields,  and  obtain  a  converged  solution  with 
a  k-e  turbulence  model.  The  'equilibrium  turbulence  model'  assumes  that 
the  production  and  dissipation  of  turbulence  kinetic  energy  are  balanced 
and  that  the  turbulent  viscosity,  p?,  can  be  calculated  from  the  mixing 
length  model,  i.e.,  Eq.  (27).  Thus  by  setting  the  source  term  of  Eq.  (27) 
to  zero  and  by  using  Eqs.  (23)  and  (25),  'equilibrium'  values  of  k  and  e 
can  be  obtained,  viz., 
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In  this  study  the  Joncs-Launder  (Ref.  16)  formulation  for  is  used 
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where  Re^  is  the  turbulence  Reynolds  number  defined  by 


(64) 


Therefore 


(65) 


Eq.  (65)  is  -a  transcendental  equation  in  Cy  (as  an  initial  guess  uj  is 
calculated  from  the  mixing  length  model  and  the  laminar  viscosity  uj  Is 
known)  which  is  solved  by  a  straightforward  ’Newton-Raphson*  iteration 
technique.  Once  this  is  done  the  value  of  Cy  can  be  substituted  into 
Eqs.  (61)  and  (62)  and  initial  guesses  of  the  k  and  e  fields  obtained. 
Boundary  conditions  for  the  k  and  e  equations  are  as  follows:  (1)  on 
inlets  the  values  of  k  and  e  are  frozen  at  their  equilibrium  values  from 
the  converged  mixing  length  solution,  (2)  on  walls  the  values  of  k  and  e 
are  set  to  aero,  (3)  on  axis  or  plane  of  symmetry  the  normal  derivatives  of 
k  and  e  are  set  to  zero,  and  (4)  on  the  exit  plane  the  second  derivatives 
of  k  and  e  are  set  to  zero. 


2.4  Numerical  Procedure 

The  numerical  procedure  used  to  solve  the  governing  equations  is  a 
consistently  split  linearized  block  implicit  (LBI)  scheme  originally 
developed  by  Briley  and  McDonald  (Ref.  12).  A  conceptually  similar  scheme 
has  been  developed  for  two-dimensional  MHD  problems  by  Llndemuth 
and  Killeen  (Ref.  24).  More  recently  Beam  and  Warming  (Ref.  25)  have 
derived  this  and  other  related  schemes  by  the  method  of  approximate 
factorization.  The  procedure  is  discussed  in  detail  in  Refs.  12  and  26. 

The  method  can  be  briefly  outlined  as  follows:  the  governing  equations  are 
replaced  by  an  implicit  time  difference  approximation.  Terms  involving 
nonllnearltles  at  the  implicit  time  level  are  linearized  by  Taylor 
expansion  in  time  about  the  solution  at  the  known  time  level,  and  spatial 
finite  difference  approximations  are  introduced.  The  result  is  a  system 
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of  multidimensional  coupled  (but  linear)  difference  equations  for  the 
dependent  variables  at  the  unknown  or  implicit  time  level.  To  solve  these 
difference  equations,  the  Douglas-Gunn  (Ref.  27)  procedure  for  generating 
alternating-direction  implicit  (ADI)  schemes  as  perturbations  of 
fundamental  implicit  difference  schemes  is  introduced  in  its  natural 
extension  to  systems  of  partial  differential  equations.  This  technique 
leads  to  systems  of  coupled  linear  difference  equations  having  narrow 
block-banded  matrix  structures  which  can  be  solved  efficiently  by  standard 
block-elimination  methods. 

The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for-the  integration  of  initial-value  problems.  The  linearization 
technique,  which  requires  an  Implicit  solution  procedure,  permits  the 
solution  of  coupled  nonlinear  equations  in  one  space  dimension  (to  the 
requisite  degree  of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no 
iteration  is  required  to  compute  the  solution  for  a  single  time  step,  and 
since  only  moderate  effort  is  required  for  solution  of  the  implicit 
difference  equations,  the  method  is  computationally  efficient.  This 
efficiency  is  retained  for  multidimensional  problems  by  using  what  might  be 
termed  block  ADI  techniques.  The  method  is  also  economical  in  terms  of 
computer  storage,  in  its  present  form  requiring  only  two  time-levels  of 
storage  for  each  dependent  variable.  Furthermore,  the  block  ADI  technique 
reduces  multidimensional  problems  to  sequences  of  calculations  which  are 
one  dimensional  in  the  sense  that  easily-solved  narrow  block-banded 
matrices  associated  with  one-dimensional  rows  of  grid  points  are  produced. 

A  more  detailed  discussion  of  the  solution  procedure  as  discussed  by 
Briley,  Buggeln  and  McDonald  (Ref.  28)  is  given  in  the  Appendix  B. 

2.5  Artificial  Dissipation 

One  major  problem  to  be  overcome  in  calculating  high  Reynolds  number 
flows  using  the  Navier-Stokes  equations  Is  the  appearance  of  spatial 
oscillations  associated  with  the  so-called  central  difference  problem. 

When  spatial  derivatives  are  represented  by  central  differences,  high 
Reynolds  number  flows  can  exhibit  a  eaw  tooth  type  oscillation  unless  some 
mechanism  is  added  to  the  equations  to  supptess  their  appearance.  This 
dissipation  mechanism  can  be  added  implicitly  to  the  equations  via  the 
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spatial  difference  molecule  (e.g.  one-sided  differencing)  or  explicitly 
through  addition  of  a  specific  term*  The  present  authors  favor  this  latter 
approach  for  two  reasons.  First,  if  a  specific  artificial  dissipation  term 
is  added  to  the  equations,  it  is  clear  precisely  what  approximation  is 
being  made.  Secondly,  if  a  specific  term  is  added  to  suppress 
oscillations,  the  amount  of  artificial  dissipation  added  to  the  equations 
can  be  easily  controlled  in  magnitude  and  location  so  as  to  add  the  minimum 
amount  necessary  to  suppress  spatial  oscillations.  Studies  can  also  be 
easily  performed  to  evaluate  the  effect  of  the  explicitly  added  dissipation 
on  the  solution. 

Various  methods  of  adding  artificial  dissipation  were  Investigated  in 
Ref.  12,  and  these  were  evaluated  in  the  context  of  a  one-dimensional  modiel 
problem.  The  model  problem  used  was  one-dimensional  flow  with  heat 
transfer.  Flow  was  subsonic  at  the  upstream  boundary,  accelerated  via  heat 
sources  until  a  Mach  number  of  unity  was  reached  and  then  accelerated  to 
supersonic  velocity  by  heat  sinks.  The  exit  back  pressure  was  raised  to 
cause  a  shock  to  appear  in  the  .supersonic  region.  This  basic 
one-dimensional  problem  contained  many  relevant  features  including  strong 
accelerations  and  the  appearance  of  a  normal  shock  wave.  Therefore,  it 
served  as  a  good  test  case  for  various  forms  of  artificial  dissipation 
which  could  be  used  in  the  presence  of  shock  waves. 

The  results  of  the  Ref.  29  investigation  led  to  the  conclusion  that 
for  the  model  problem  a  second  order  artificial  dissipation  approach  was 
the  best  of  those  considered.  This  approach  adds  a  term  of  the  form 
yart  32$/322  a/3Z  {vart  34/3Z}  to  each  governing  equation  where 
♦  ■  p,  u,  v,  w  for  the  continuity,  x-momentum,  y-raomentum  and  z-momentum 
equations  respectively  and  vart  is  determined  by 


UjAZ 


(66) 


The  A 7,  in  Eq.  (66)  is  the  distance  between  grid  points  In  n  given 
coordinate  direction,  Ug  is  the  velocity  in  this  direction,  og  is  the 
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artificial  dissipation  parameter  for  this  direction  and  v  is  the  effective 
kinematic  viscosity.  The  equation  determines  vart  with  vart  taken  as 
the  smallest  non-negative  value  which  will  satisfy  the  expression.  It 
should  be  noted  that  in  two  space  dimensions  each  equation  contains  two 
artificial  dissipation  terns,  one  in  each  coordinate  direction.  For 
example,  the  6treamwise  momentum  equation  expressed  in  two-dimensional 
Cartesian  coordinates  would  contain  the  artificial  dissipation  terms 
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3.0  RESULTS 


The  primary  objective  of  the  computational  effort  is  to  demonstrate 
the  capability  of  the  previously  described  numerical  technique  to 
accurately  calculate  the  flow  in  realistically  configured  labyrinth  seals 
at  typical  operating  conditions*  It  is  hoped  that  the  achievement  of  this 
objective  will  provide  insight  into  the  details  of  the  flow  structure  not 
easily  obtainable  by  experimental  methods,  and  that  this  demonstration  will 
lead  to  the  use  of  this  analytical  tool  in  the  design  of  future  labyrinth 
seals*  In  order  to  achieve  the  above  objective,  a  series  of  calculations 
were  performed.  These  calculations  can  logically  be  divided  into  four 
categories:  (1)  calculations  designed  to  show  the  ability  of  the  numerical 
procedure  to  accurately  predict  the  leakage  of  various  labyrinth  seal 
designs  at  typical  operating  conditions,  (2)  calculations  to  obtain  a 
performance  curve  for  two  seal  configurations,  (3)  calculations  of  the  flow 
for  two  seal  configurations  which  can  be  compared  with  experimental  data 
obtained  at  Allison  Gas  Turbine  Operations  (Ref.  30)  for  the  same  seal 
configurations  and  at  the  same  operating  conditions,  and  (A)  a 
demonstration  calculation  of  the  ability  of  the  numerical  procedure  to 
calculate  the  flow  in  a  sample  configuration  while  the  seal  assembly  Is 
rotating. 

Before  going  Into  the  details  of  the  calculations.  It  Is  desirable  to 
define  the  labyrinth  seal  nomenclature  thut  Is  appropriate  to  the 
configurations  considered  In  this  study.  The  geometric  capability  of  the 
computer  program  developed  under  this  study  allows  for  the  analysis  of  a 
wide  variety  of  both  conventional  straight-through  and  stepped  seal 
configurations.  Nomenclature  that  is  common  to  both  types  of  seals  Is 
shown  In  Pig*  (2).  CL  Is  the  clearance  between  the  seal  knife  tlp(s)  and 
the  stator  or  land.  KT  represents  the  knife  blade  thickness  at  the  dp. 

KB  ts  the  knife  blade  taper  angle,  and  K6  Is  the  knife  blade  slant  angle 
relative  to  the  rotor  surface.  Kll  represents  the  knife  blade  height  and 
KP  ts  the  distance  between  successive  knife  blades.  For  stepped  seal 
configurations  DTC  is  used  to  represent  the  minimum  horizontal 
dlstanee-to-contact  between  the  knife  blade  and  the  stator,  and  SH  la  the 
step  height.  In  addition  for  stepped  seal  configurations  the  leakage  flow 
direction  is  referred  to  as  either  from  the  large-to-small  seal  diameter 
(LTSD)  or  trom  the  small-to-large  seal  diameter  (STL0). 
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The  flow  reference  conditions  were  determined  in  a  similar  manner  for 
all  cases  considered.  The  inlet  stagnation  pressure,  P0,  and  stagnation 
temperature,  T0,  were  specified  as  was  the  ratio  of  upstream  tc 
downstream  pressure,  P0/Pd*  In  addition  an  estimate  of  the  leakage 
rate  (mass  flow),  W,  through  the  seal  was  obtained  either  from  the  results 
of  the  Allison  experimental  program  or  the  Allison  design  model.  Given  the 
upstream  pressure  and  temperature,  the  upstream  density,  p0,  can  be 
calculated  from  the  perfect  gas  law  while  the  upstream  viscosity,  y,  can  be 
calculated  from  Sutherland's  viscosity  law,  Ref.  31.  The  average  inlet 
velocity,  V,  can  be  calculated  from  the  relationship  (for  the  cases 
considered  the  inlet  Mach  numbers  were  on  the  order  of  0.01,  hence  the 
stagnation  and  the  static  conditions  are  virtually  identical) 


where  A  is  the  known  inlet  area.  All  cases  considered  in  this  study  had  a 
spanwise  distance,  of  6.28  in.  ■  .160  m.  The  above  now  yields  enough 
Information  so  that  reference  (or  upstream)  Reynolds  and  Mach  numbers  can 
be  calculated.  Using  the  above  upstream  variables  as  reference  quantities, 
the  governing  equations  can  be  nondlmenslonallzed.  When  a  case  is 
converged  for  a  given  ratio  of  upstream  to  downstream  pressure  ratio 
(sometimes  referred  to  as  the  seal  expansion  ratio  r  -  P0/PdK  the 
leakage  rate  through  the  seal  can  then  be  recalculated  and  compared  to  the 
Initial  estimate.  For  instance  if  the  calculated  nondlmenslonal  average 
upstream  velocity  V  is  0.5,  the  computed  leakage  rate  would  be  50%  of  the 
estimated  rate.  The  calculated  reference  Reynolds  and  Mach  numbers  would 
correspondingly  have  to  be  reduced  by  50%.  The  Initial  conditions  for  all 
cases  considered  was  to  assume  that  the  flow  was  Initially  at  rest  at  the 
stagnation  conditions.  The  back  pressure  was  then  gradually  lowered  until 
the  desired  value  of  Fp  was  obtained.  The  solution  was  then  time  marched 
until  a  steady  state  solution  was  obtained*  The  output  of  each  calculation 
consists  of  the  fields  of  the  independent  and  dependent  variables,  l.c., 
the  Cartesian  or  cylindrical  polar  coordinates  and  the  velocity  components, 
density,  enthalpy  and.  If  appropriate,  the  turbulence  kinetic  energy  and 
the  turbulence  dissipation,  tn  addition  any  of  the  derived  variables  such 
as  pressure,  temperature,  viscosity,  stream  function,  total  temperature  and 
pressure,  etc.  can  be  calculated  and  displayed. 
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3.1  Calculation  of  Leakage  Rates 


The  majority  of  the  calculations  performed  under  this  effort  were 
done  in  the  first  category  of  calculations  and  were  designed  to  demonstrate 
the  ability  of  the  computational  procedure  to  predict  the  leakage  rate  for 
a  large  number  of  straight-through  and  stepped  seal  configurations  at 
various  flow  conditions.  A  total  of  18  cases  were  run,  13  straight-through 
seal  configurations  and  5  stepped  seal  configurations.  Tables  1  and  2  and 
Figs.  3  and  4  give  a  brief  synopsis  of  the  configurations  considered,  the 
flow  conditions  and  the  modelling  assumptions  under  which  the  calculations 
were  performed.  Several  of  the  calculations  were  made  before  the  energy 
equation  option  was  operational  in  the  computer  code.  For  these 
calculations  the  energy  equation  was  approximated  by  assuming  that  the 
stagnation  enthalpy,  h0,  was  constant.  For  flow  in  the  Mach  number  range 
considered  in  this  study  (essentially  incompressible  on  the  inlet  plane  to 
a  peak  in  the  flow  field  of  between  0.7  and  l.b  depending  on  the  flow 
conditions  and  the  geometry  of  the  seal),  this  assumption,  which  neglects 
heat  conduction  and  stress  work,  is  a  reasonable  approximation  to  the 
physics*  If,  however,  the  walls  are  nonadlabatlc,  i.e.,  either  highly 
cooled  or  heated,  the  assumption  of  zero  heat  transfer  is  Invalid.  For  the 
cases  considered  under  this  study,  the  walls  were  neither  heated  nor 
cooled.  Hence  the  adiabatic  assumption  was  valid.  In  addition  the  earlier 
calculations  presumed  that  the  inlet  flow  was  fully  developed,  and,  hence. 
In  these  calculations  a  one-seventh  (1/7)  power  law  (Ref.  22)  was  used  for 
the  inlet  streamwlse  velocity  profile.  For  later  calculations,  a  boundary 
layer  equal  to  SQZ  of  the  clearance  was  assumed  on  both  the  Inlet  section 
of  the  rotor  and  the  land.  Within  the  boundary  layers,  the  method  of  Haisc 
and  McDonald  (Ref.  23)  was  used  to  obtain  the  turbulent  streamwlse  velocity 
profile.  The  mixing  length  turbulence  model  utilised  was  a  hybrid  model 
consisting  of  a  Williamson's  model  (Ref.  32)  in  the  regions  away  from  solid 
wails  and  a  van  Driest  damped  model  (Ref.  33)  in  the  regions  near  walls. 

In  practice,  both  mixing  length  values  were  calculated  at  each  grid  p>int 
and  the  locally  minimum  value  chosen.  This  gives  a  smooth  variation  of 
mixing  length  throughout  the  flow  field.  In  the  later  calculations,  the 
flows  were  calculated  with  both  the  mixing  length  and  the  previously 
discussed  k-c  turbulence  model. 
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The  number  of  grid  points  utilized  for  each  calculation  is  presented 
in  Tables  1  and  2.  The  basic  philosophy  was  to  concentrate  the  grid  points 
in  the  regions  where  the  largest  physical  gradients  of  the  dependent 
variables  were  expected.  These  areas  were:  (1)  in  the  region  between  the 
tip  of  the  knife  blades  and  the  land,  (2)  in  the  wall  boundary  layers, 
and  (3)  in  the  vicinity  of  rapid  expansion  or  compression.  The  total 
streamwise  extent  of  the  physical  domain  chosen  for  the  computations  was 
300  clearances  in  length.  The  inlet  plane  was  chosen  at  a  distance  of 
50  clearances  upstream  of  the  front  face  root  of  the  first  knife  leaving  as 
much  as  250  clearances  (depending  on  the  seal  configuration)  from  the  last 
knife  to  the  exit  plane.  The  relatively  large  extent  of  the  domain 
downstream  of  the  last  knife  was  required  by  the  existence  of  a  large 
streamwise  recirculation  zone  downstream  of  the  last  knife  for  each 
case  investigated.  This  large  domain  was  needed  to  ensure  that  the 
recirculation  zone  would  remain  inside  the  region  chosen  for  the 
computation.  A  typical  distribution  of  grid  points  is  shown  in  the 
vicinity  of  the  seal  assembly  (see  Fig.  5).  In  the  regions  both  upstream 
and  downstream  of  the  knives,  the  streamwise  grid  spacing  was  considerably 
larger  since  in  these  regions  the  streamwise  gradients  were  relatively 
small. 

Both  Tables  1  and  2  compare  the  calculated  leakage  rates,  Wcaic, 
with  the  Allison  correlation  leakage  rates,  Wcorr.  The  Allison 
correlation  leakage  rates  are  a  composite  of  the  calculated  leakage  rates 
as  determined  by  the  Allison  design  model  and  experimentally  measured 
leakage  rates.  The  Allison  design  model  (Ref.  30)  is  based  on  a  multiple 
regression  analysis  of  a  large  bank  of  experimental  data.  The  design  model 
can  predict  the  leakage  rates  for  a  wide  variety  of  seal  geometries  at 
various  flow  conditions.  The  data  used  as  the  basis  for  this  model  are 
taken  from  experiments  performed  on  full-scale  labyrinth  seals  since  this 
model  is  Intended  to  be  used  to  design  contemporary  gas  turbine  engine 
seals.  The  calculations  performed  under  this  effort,  however,  were  made 
for  large-scale  (nominally  10X  scale  for  the  straight-through  seals  and 
5X  scale  for  the  steppod  seals)  labyrinth  seals  as  were  the  leakage  rate 
measurements  performed  under  the  experimental  portion  of  this  effort.  When 
an  experiment  had  been  performed  on  the  particular  labyrinth  seal  at  the 
specified  flow  conditions,  the  Allison  correlation  leakage  rate  was  taken 
to  oe  the  experimentally  determined  value.  This  is  noted  on  Tables  1 
and  2  by  the  Allison  designated  test  number.  However  for  several  of  the 
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cases  calculated,  no  experimental  leakage  rates  were  available.  To 
determine  correlation  leakage  rates  for  these  cases,  the  following 
procedure  was  used.  First  the  cases  were  divided  into  straight  through  and 
stepped  seal  configurations.  Tables  1  and  2,  respectively*  For  the 
straight  through  seals,  the  Allison  design  model  was  used  to  predict  the 
leakage  rates  for  the  cases  for  which  experimental  data  existed,  i.e.  cases 
5,  6,  7  and  7A.  For  these  cases  it  was  found  that  the  Allison  design  model 
underpredicted  the  leakage  rates  by  an  average  of  12Z.  It  is  felt  that 
this  underprediction  is  due  to  the  Reynolds  number  effect  of  the  smaller 
full-scale  seals  used  in  the  design  model.  Hence,  to  account  for  this 
effect  for  the  cases  where  experimental  data  were  not  available,  the 
Allison  design  model  was  used  to  determine  preliminary  values  which  were 
then  scaled  by  the  factor  of  I « 12  to  account  for  the  large-scale  seals  used 
in  the  calculation. 

Of  the  four  stepped  seal  cases  considereu  in  this  study,  experimental 
data  existed  for  two  cases,  cases  12  and  13.  These  values  were  used  as  the 
Allison  correlation  leakage  rates,  Wca}c,  for  these  cases.  For  case  11 
(a  lOX-scale  model),  experimental  results  existed  for  a  3X-scale  seal  at 
the  same  flow  conditions.  This  result  was  compared  with  the  Allison  design 
model  prediction  for  the  full-scale  seal.  The  ratio  of  leakage  rates, 

1.09,  was  applied  to  the  lOX-scale  Allison  design  model  leakage  rate  to 
obtain  the  final  correlation  value*  For  case  14,  the  experimental  leakage 
rate  for  case  12  (the  three  knife  stepped  seal  with  rectangular  knives, 
STLO)  was  compared  with  the  design  model  full-scale  Leakage  rate  for  the 
case  12  seal.  The  ratio  of  these  two  leakage  rates,  0.96,  was  then  used  to 
correct  the  design  model  leakage  rate  for  case  14.  The  assumption  was 
that  the  ratio  between  the  experimental  Leakage  and  the  design  model  for 
the  tapered  knives  is  the  same  as  that  for  rectangular  knives.  The  results 
for  the  stepped  seal  cases  arc  contained  in  Table  2. 

Examination  of  the  results  for  the  straight-through  seal  cases 
Table  1  shows  generally  good  agreement  between  the  calculated  leakage 
rates  and  Allison  correlation  leakage  rates.  To  better  understand  the 
implications  of  the  results  presented  in  Table  1,  it  is  perhaps  desirable 
to  examine  the  cases  Individually  and/or  in  logical  groupings. 

Cases  for  which  cxperl&enal  data  arc  available,  excepting  case  8,  the 
three-knife  straight-through  seal  with  tapered  knives,  show  excellent 
agreement.  The  average  discrepancy  is  3Z.  The  calculated  leakage  rate 
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for  case  8  is  from  15  to  25  percent  lower  than  that  measured  depending  on 
whether  a  mixing  length  or  a  k-e  turbulence  model  was  utilized  in  the 
calculation..  For  this  case  the  k-e  turbulence  model  gives  larger  values  of 
turbulent  viscosity  than  the  mixing  length  model  and,  hence,  thicker 
boundary  layers  and  less  leakage.  On  the  other  hand  the  predicted  and 
measured  leakage  rates  for  case  7,  the  three  knife  straight-through  seal 
with  tapered  knives  at  a  pressure  ratio  of  4.98,  shows  excellent  agreement, 
the  calculated  value  predicting  only  1.5  percent  more  than  the  measured 
value.  The  flow  for  this  case  was  choked  as  the  peak  Mach  number  was 
1.398.  Examination  of  the  flow  for  the  other  straight-through  case  where 
the  flow  was  choked  and  experimental  data  existed,  case  5,  indicates  a 
general  ability  to  accurately  calculate  leakage  rates  for  choked  cases; 
the  discrepancy  between  the  calculated  and  measured  leakage  rates  for  this 
case  is  approximately  0.9  percent. 

For  case  1,  the  single  knife  straight-through  seal  with  a  rectangular 
knife  at  a  pressure  ratio  of  3.98,  no  experimentally  measured  value  of 
leakage  rate  was  available.  The  correlation  value  predicts  approximately 
30  percent  more  leakage  than  the  calculated  value.  This  might  be  due  to 
the  use  of  only  21  grid  points  in  the  transverse  direction  which  could 
possibly  lead  to  an  underprediction  of  mass  flux  (this  was  the  first  case 
run  under  this  effort). 

The  three-knife  straight-through  seal  with  rectangular  knives  run  at  a 
pressure  ratio  of  2.0,  case  2,  shows  a  reasonable  agreement  between  the 
calculated  and  the  correlation  leakage  rates  (no  experimental  data  were 
available  for  this  case).  The  calculated  value  of  leakage  was  1.7  percent 
lower  than  the  correlation. 

For  the  cases  of  the  single  knife  and  triple  knife  straight-through 
seals  with  slanted  knives  at  a  pressure  ratio  of  2.0,  (cases  3  and  4)  the 
agreement  between  the  calculated  leakage  rates  and  the  correlation  values 
are  reasonable  with  an  ovorpredlction  of  leakage  rates  of  12  and  7  percent 
for  cases  3  and  4  respectively.  From  a  computational  viewpoint,  these 
cases  are  extremely  difficult  because  of  the  geometric  requirements  that 
the  width  of  the  base  of  the  knife  blade  Is  large  compared  to  the  width 
of  the  tip  of  the  knife  blade,  i.e.,  a  ratio  on  the  order  of  100:1. 

This  not  only  causes  difficulty  in  generating  a  grid  structure,  but 
inherently  yields  a  coordinate  system  that  is  highly  skewed  (as  opposed 
to  orthogonal)  relative  to  the  other  cases  considered  in  this  study. 

Because  the  flow  running  up  the  leading  edge  surface  of  the  knife 


is  oriented  at  an  angle  greater  than  90°  to  the  incoming  flow,  this  results 
in  extremely  large  gradients  in  the  vicinity  of  the  leading  edge  of  the  tip 
of  the  knife,  i.e.,  in  the  entrance  region  of  the  gap  between  the  first 
knife  and  the  land.  The  large  gradients  require  the  use  of  a  large  number 
of  grid  points  in  the  vicinity  of  the  leading  edge  in  order  to  resolve  the 
large  gradients.  The  existence  of  the  large  gradients  in  this  region  leads 
to  strong  dissipation  and  hence  the  excellent  performance  of  this  seal. 

In  the  case  of  the  worn  single  knife  straight-through  seal,  the 
calculated  leakage  rate  underpredicts  the  measured  rate  by  5.8  percent  and 
1.2  percent  for  the  mixing  length  and  a  k-e  turbulence  modelsrespectively. 
In  this  case  the  k-e  turbulence  model  in  general  predicts  slightly  lower 
values  of  turbulent  viscosity  and  hence  a  slightly  higher  leakage  rate.  In 
both  cases  the  discrepancy  between  the  measured  and  calculated  rates  are 
low. 

The  calculation  of  the  flow  in  the  three  knife  straight-through  seal 
with  tapered  knives  at  a  pressure  ratio  of  2.0,  case  8,  was  run  in  three 
modes:  (1)  a  mixing  length  model  with  a  fully  developed  1/7  power  law 
turbulent  inlet  profile  with  a  51  x  71  computational  grid,  (2)  a  mixing 
length  model  with  a  turbulent  inlet  profile  on  the  inlet  land  and  rotor 
equal  to  50%  of  a  clearance  height  with  a  61  x  101  grid  and  (3)  a  k-e 
turbulence  model  with  the  same  inlet  profile  and  grid  structure  as  in  (2). 
For  the  first  two  modes,  the  calculated  leakage  rate  was  approximately  15 
percent  less  than  the  measured  value,  while  for  the  k-e  turbulence  model 
the  calculated  leakage  rate  was  approximately  25  percent  lower.  In  light 
of  the  excellent  agreement  for  the  three  knife  straight-through  seal  with 
rectangular  knives  and  the  single  knife  straight  through  seal  with  a 
tapered  knife,  where  the  physical  processes  should  be  similar  to  the  three 
knife  straight-through  seal  with  taoered  knives,  the  amount  of  the 
discrepancy  was  unexpected.  Possible  reasons  for  the  discrepancy  could  be 
the  existence  of  leaks  in  the  experimental  apparatus,  three-dimensional 
effects  in  the  apparatus,  the  sudden  expansion  downstream  of  the  last  knife 
blade  (which  did  not  exist  in  the  calculation)  or  numerical  truncation 
error  in  the  calculation.  Without  either  repeating  the  experiment  and/or 
the  calculation,  it  is  difficult  to  draw  any  firm  conclusions  for  case  8. 

In  addition,  it  is  interesting  to  note  that  pressure  taps  were  placed  at 
key  locations  In  the  experimental  apparatus  (see  Fig.  6),  and  that  the 
agreement  between  the  pressure  calculation  and  experiment  were  excellent. 
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For  the  three  knife  straight-through  seal  with  tapered  knives  and 
with  a  rough  land  (the  roughness  was  produced  by  attaching  a  30  grit 
sand  paper  to  the  land)  at  a  pressure  ratio  also  of  2.0  (case  9), 
the  agreement  between  the  calculated  and  measured  leakage  rate  is 
reasonable.  The  calculated  leakage  value  is  1.4  percent  lower  than  the 
measured  value.  The  wall  roughness  was  simulated  by  using  the  method  of 
van  Driest  (Ref.  33)  where  a  slip  velocity  is  assumed  on  the  rough  wall. 
Computationally  this  has  the  effect  of  making  the  boundary  layer  thinner 
and  hence  more  mass  can  pass  through  the  device.  The  initial  condition  for 
this  case  was  taken  as  the  converged  solution  for  case  8B«  The  no-slip 
condition  on  the  land  was  then  replaced  by  the  van  Driest  model,  and  a 
converged  solution  obtained.  It  is  interesting  to  note  that  the  effect  of 
the  roughness  was  to  increase  the  peak  Mach  number  from  0.776  to  0.946. 

Finally  the  last  case  considered  was  the  three  knife  straight-through 
seal  with  tapered  knives  and  with  injection  at  a  pressure  of  2.0 
(case  10).  The  injection  rate  was  chosen  as  nominally  10  percent  of  the 
total  leakage  rate.  The  injection  port  was  positioned  in  the  land  midway 
between  the  first  and  second  knives.  Again  the  initial  condition  was  taken 
as  the  converged  result  of  case  8b.  The  injection  had  the  same  effect  as 
land  roughness  on  the  overall  leakage  rate.  The  effect  of  the  injection 
(from  a  computational  viewpoint)  is  that  the  amount  of  leakage  is  increased 
over  the  noninjection  case.  (Note  that  no  correlation  value  of  leakage 
rate  is  available  for  this  case).  One  possible  explanation  is  that  the 
injected  fluid  la  not  exposed  to  the  large  losses  associated  with  the 
leading  edge  of  the  first  knife.  The  injection  has  only  the  minimal  local 
effect  of  forcing  the  flow  emerging  from  between  the  first  knife  and  land 
slightly  further  into  the  cavity  region  between  the  first  and  second 
knives.  Thus  the  tentative  conclusion  is  that  for  this  injection  case  is 
the  effect  is  negative,  i.e.,  the  leakage  rate  is  Increased.  Perhaps  the 
more  desirable  location  for  Injection  would  be  in  the  vicinity  of  the 
leading  edge  of  the  first  knife  to  further  increase  the  losses  in  this 
critical  region. 

The  results  for  the  stepped  seal  configurations  are  presented  In 
Tahle  2.  All  canes  were  run  In  the  flow  direction  of  saall-to-large  seal 
diameters  (STM)  mode).  Kxpertmentai  data  existed  for  cases  12  and  13. 

Cases  II  and  14  had  the  predicted  correlation  leakage  rates  previously 


discussed.  For  the  stepped  seal  configurations,  only  case  11  was  run 
under  choked  conditions.  For  this  case  the  pressure  ratio  was  2.5 
(all  other  stepped  configurations  were  run  at  a  pressure  ratio  of  2.0), 
and  the  peak  calculated  Mach  number  was  1.229.  The  calculated  leakage 
rate  was  approximately  12  percent  larger  than  the  correlated  value. 

A  similar  trend  was  also  noted  for  case  12,  the  three  knife  stepped  seal 
with  rectangular  knives,  where  the  predicted  leakage  rate  was  25  percent 
higher  than  the  experimentally  determined  rate.  Possible  reasons  for  this 
discrepancy  would  be  as  was  previously  discussed  for  the  three  knife 
straight-through  seal  with  tapered  knives.  Again  it  is  of  interest  to  note 
that  the  correlation  between  the  measured  and  calculated  static  pressures 
at  key  points  in  the  seal  Is  excellent  (see  Fig.  7).  This  result  is 
similar  to  that  observed  for  the  case  of  the  three  knife  straight  through 
seal  with  tapered  knives.  The  other  stepped  seal  configuration  for  which 
experimental  data  existed  was  case  13,  the  single  knife  stepped  seal  with  a 
slanted  knife  configuration.  In  this  case,  the  predicted  leakage  rate  was 
8  percent  higher  than  the  measured  value.  The  three  knife  stepped  seal 
with  tapered  knives  configuration  (case  14)  was  calculated  with  both  a 
mixing  length  and  a  k-s  turbulence  model.  The  calculated  results  were  not 
significantly  different.  The  mixing  length  model  predicting  a  leakage  rate 
about  10  percent  higher  than  the  correlated  value,  and  the  k-e  turbulence 
model  predicting  a  leakage  rate  about  11  percent  higher. 

In  general  the  capability  of  the  numerical  procedure  to  accurately 
predict  the  leakage  rates  for  a  wide  variety  of  seal  configurations  under 
various  flow  conditions  appears  to  be  juslfled.  In  the  worst  case  the 
prediction  differed  from  the  Allison  correlation  value  by  30  percent.  In 
the  vast  majority  of  cases  a  more  typical  variation  would  be  5-10  percent, 
and  in  many  cases  the  predicted  and  measured  values  differed  by  only  a  few 
percent.  Since  this  was  the  first  effort  in  applying  this  numerical 
procedure  to  labyrinth  seals,  it  Is  to  be  expected  that  future  results 
would  be  even  more  encouraging.  A  primary  advantage  of  the  calculation 
procedure  is  that  details  of  the  flow  are  calculated  at  every  point  In  the 
flow,  hence  the  manner  in  which  the  flow  develops,  the  basic  flow  patterns, 
etc.  can  be  discerned  immediately.  Typical  CPU  calculation  times  for  these 
cases  were  on  the  order  of  1.5  hours  on  a  CDC  7600  at  the  Ballistic 
Research  Laboratory,  l.e*,  about  $500  per  case  at  the  overnight  rate.  In 


37 


addition,  recent  efficiencies  incorporated  into  the  computer  program  used 
in  the  above  calculations,  viz.,  in  the  vectorization  on  the  CRAY 
computers,  have  resulted  and  will  result  in  even  lower  cost  per  case  rates 
in  the  future. 

Plots  of  streamlines,  velocity  vectors  and  Mach  number  contours  are 
shown  for  representative  cases  calculated  in  the  leakage  rate  study 
(Figs.  8  to  13).  When  the  earlier  cases  were  run,  the  capability  to 
calculate  streamlines  was  not  operational.  However,  basic  streamline 
patterns  can  be  inferred  from  the  velocity  vector  plots.  The  plots 
represent  the  flow  only  in  the  region  of  the  knives.  For  all  cases 
investigated,  a  large  clockwise  streamwlse  recirculation  was  generated 
downstream  of  the  last  knife  for  flow  from  left  to  right.  Usually  this  v?as 
accompanied  by  a  small  counterclockwise  recirculation  zone  at  the  junction 
of  the  downstream  base  of  the  last  knife  at  the  rotor.  In  some  cases  small 
clockwise  recirculation  zones  were  calculated  in  the  gap  between  the 
leading  edge  of  the  first  knife  and  the  land.  When  no  recirculation  zone 
was  calculated  in  this  region,  the  flow  was  significantly  decelerated  due 
to  the  strong  adverse  pressure  gradient  in  this  region.  In  all  cases  a 
recirculation  pattern  existed  In  the  cavity  regions.  In  some  cases  a  snail 
counterclockwise  recirculation  zone  would  exist  at  the  trailing  edge  of  the 
knives  due  to  the  separation  of  the  flow.  At  the  junction  of  the  upstream 
face  of  the  first  knife  and  the  rotor  cylinder,  a  small  clockwise 
recirculation  zone  would  usually  exist*  In  addition,  in  the  case  of 
stepped  seals,  the  flow  would  separate  off  the  convex  corner  of  the  land 
step  forming  a  counterclockwise  recirculation  pattern.  In  general,  the 
controling  location  of  the  flow  for  the  cases  investigated  appeared  to  be 
the  entrance  region  between  the  tip  of  the  first  knife  and  the  land.  In 
this  region  the  streamwlse  flow  would  accelerate  and  large  (on  the  order  of 
the  streamwlse  velocity)  transverse  flows  would  be  generated.  At  the  inlet 
regions  of  any  downstream  knives,  the  accc  orations  (and  hence  lossea) 
would  be  much  loss  as  the  flow  did  not  significantly  turn  into  the  cavity 
regions.  Thus,  the  flow  enters  the  clearance  gaps  of  th<  subsequent  knives 
with  relatively  small  tansverse  velocities  and,  hence,  the  loss  is  not 
nearly  as  great  as  at  the  first;  knife.  For  the  cases  where  the  flow  is 
choked,  an  expansion  pattern  was  predicted  in  the  region  of  the  exit  of  the 
last  knife.  This  can  be  seeu  in  the  hach  number  contour  plots  where  the 
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gradient  in  Mach  number  is  large,  i.e.,  the  Mach  number  contours  are  close 
together. 

3.2  Generation  of  Performance  Curves 


The  purpose  of  the  second  classification  of  calculation  was  to 
demonstrate  that  the  computer  code  could  be  used  to  generate  designers’ 
curves  or  performance  plots  for  specific  seal  geometries.  Specifically,  it 
is  desired  to  construct  curves  of  the  flow  parameter,  $, 


<*>  = 


(69) 


(where  again  W  is  the  leakage  rate,  T0  is  the  upstream  stagnation 
temperature,  P0  is  the  upstream  stagnation  pressure  and  A  is  the 
clearance  area  over  one  of  the  knives)  versus  the  expansion  ratio, 

Pq/Pd,  for  a  specific  seal  geometry.  In  this  study  two  seal  geometries 
were  considered:  (1)  a  three  knife  straight  through  seal  with  tapered 
knives.  Pig.  3f,  and  (2)  a  three  knife  stepped  seal  with  tapered  knives. 
Fig.  4c.  Both  of  these  geometries  were  considered  in  the  previously 
described  leakage  calculations.  The  main  idea  here  is  to  generate  the 
curves  for  cases  for  which  experimental  data  exist  and  for  which  the 
Allison  design  model  can  be  used  to  generate  similar  curves.  Demonstration 
of  the  ability  of  the  numerical  procedure  to  produce  reasonable  performance 
plots  for  these  two  designs  would  be  part  of  the  overall  validation  process 
for  the  code.  In  addition,  it  would  lend  credence  to  the  use  of  the 
numerical  procedure  to  generate  performance  curves  for  either  advanced  seal 
configurations  for  which  data  did  not  exist  or  for  investigating  variants 
of  existing  seal  configurations  such  as  the  effects  of  injection,  the 
effects  of  various  rotational  speeds,  the  effects  of  wear  on  various  seal 
designs,  etc.  These  curves  could  then  be  used  by  a  designer  of  labyrinth 
seal  systems  to  investigate  various  candidate  configurations  without  the 
need  for  setting  up  experimental  rigs.  Using  this  as  a  method  of 
eliminating  undesirable  configurations,  the  designer  would  then  be  able  to 
experimentally  investigate  the  performance  of  the  remaining  candidate 
configurations.  Considerable  savings  could  be  realized  by  the  use  of  such 
a  process  and,  in  addition,  a  better  seal  could  be  produced. 


The  technique  utilized  in  this  study  to  obtain  the  performance  plots 
was  started  from  a  converged  solution  at  a  pressure  ratio  of  2.0.  For  the 
two  seal  configurations  considered  in  this  study  the  converged  solutions 
from  the  previously  described  leakage  investigation  (cases  8B  and  14A)  were 
used  as  the  initial  conditions,  i.e.t  the  pressure  ratio  of  2.0  cases.  The 
back  pressures  for  these  two  cases  were  further  lowered  to  a  second  desired 
pressure  and  a  converged  solution  obtained  for  these  cases.  This  process 
was  continued  at  other  back  pressures.  For  this  study  converged  solutions 
were  obtained  for  pressure  ratios  of  2.0,  2.5,  3.5  and  5.0.  Initially  for 
both  seal  configurations  the  flow  was  unchoked.  However,  as  the  back 
pressure  was  lowered  the  flow  choked.  For  both  cases  considered,  a  mixing 
length  turbulence  model  was  utilized. 

Overall  the  use  of  the  computational  procedure  to  produce  a 
performance  curve  for  the  two  seal  configurations  seems  to  be  well 
justified.  The  two  curves  produced,  although  not  yielding  quantitatively 
the  same  results  as  the  design  model,  gave  curves  that  were  qualitatively 
similar  to  the  measured  data.  In  addition,  the  prediction  of  the  choke 
points  appear  to  be  In  reasonable  agreement  with  the  data.  It  is  to  be 
expected,  although  this  has  not  been  demonstrated,  that  similar  curves 
could  be  generated  for  the  other  configurations  considered  in  the  leakage 
rate  study.  In  addition,  the  performance  curves  could  also  be  generated 
for  alternate  designs  such  as  smaller  or  larger  clearance,  knife  tip 
thicknesses,  slant  angle,  etc.  A  major  advantage  would  be  the  ability  to 
generate  performance  plots  for  configurations  not  presently  in  the  Allison 
or  other  design  models. 

Performance  curves  were  generated  for  two  seal  configurations: 

(1)  the  three  knife  straight-through  seal  with  tapered  knives  and  (2)  the 
three  knife  stepped  seal  with  tapered  knives.  These  seal  configurations 
are  the  same  as  those  analyzed  for  leakage  rates  under  cases  8  and  24, 
l.e..  Fig.  3f  and  Fig*  3e,  respectively.  Results  for  these  two 
configurations  are  presented  In  Figs.  14  and  15.  For  both  cases  the 
calculated  values  of  the  flow  parameters  are  compared  with  the  experimental 
data,  the  Allison  full-scale  design  model  and  the  Allison  design  model 
corrected  for  the  effect  of  Reynolds  number  (previously  discussed).  For 
the  three  knife  stepped  seal  with  tapered  knives  no  experimental  data  were 
available,  so  In  that  case  the  experimental  data  for  the  three  knife 
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stepped  seal  with  rectangular  knives  was  used.  It  Is  expected  that  the 
tapering  of  the  knives  would  have  a  small  effect  on  the  performance  of  the 
seal. 

The  results  for  the  three  knife  straight-through  seal  with  tapered 
knives  presented  In  Fig.  14  show  that  the  calculated  values  of  flow 
parameter  are  lower  than  the  measured  values  of  flow  parameter  by 
approximately  15  per  cent  at  all  pressure  ratios.  This  is  consistent  with 
the  results  shown  in  the  leakage  study.  Possible  explanations  for  the 
difference  between  the  measured  and  calculated  flow  parameter  could  be  the 
following:  (1)  leaks  in  the  experimental  apparatus  which  would  yield  a 
larger  mass -flux,  (2)  hot  wire  measurement/calibration  errors, 

(3)  under-resolution  of  the  boundary  layers  in  the  calculation  which  would 
lead  to  thicker  boundary  layers  and  hence  less  computational  mass  flux  and, 

(4)  errors  associated  with  the  turbulence  model  used  in  the  calculation. 

In  addition,  it  is  also  possible  that  the  three-dimensional  experimental 
apparatus  might  have  had  significant  three-dimensional  effects,  and  thus 

a  discrepancy  would  exist  when  comparing  the  results  with  a  two-dimensional 
calculation.  The  shape  of  the  curves  is  similar  and  the  predicted 
choke  point  (somewhere  between  a  pressure  ratio  of  2.5  and  3.5)  is 
consistent  with  the  experimentally  observed  choke  pressure  ratio  of  2.72. 
The  full-scale  design  model  values  of  flow  parameter  are  approximately 
5  percent  higher  than  the  calculated  values  while  the  design  model  values 
corrected  for  the  apparent  effect  of  Reynolds  number  are  approximately 
15  percent  higher  than  the  calculated  values. 

The  results  for  the  three  knife  stepped  seal  with  tapered  knives  is 
presented  in  Fig.  15.  In  this  case  the  predicted  values  of  the  flow 
parameter  are  approximately  10  percent  higher  than  the  measured  values. 

The  shape  of  the  calculated  performance  curve  is  of  the  same  form  as  the 
measured  values.  The  calculated  choke  value  of  the  pressure  ratio  is 
approximately  3. 5-4.0  which  Is  consistent  with  the  measured  value  of  3.84. 
At  the  higher  values  of  pressure  ratio  there  is  some  scatter  in  the  data. 

In  addition  from  Fig.  15  it  can  be  seen  that  the  full-scale  Allison  design 
model  uniformly  overpredicts  the  flow  parameter  for  all  pressure  ratios. 


3.3  Comparison  of  Calculated  and  Experimental  Results 


To  further  validate  the  ability  of  the  computational  procedure  to 
accurately  calculate  the  flow  in  labyrinth  seals,  a  comparison  between 
available  experimental  data  and  computed  results  was  made.  Two  labyrinth 
seal  configurations  were  considered  in  this  effort:  (1)  the  three  knife 
straight  through  seal  with  tapered  knives  (see  Fig.  3f)  and  (2)  the  three 
knife  stepped  seal  with  tapered  knives  (see  Fig.  4c).  For  both  of  the 
configurations  extensive  hot  wire  measurements  were  made  as  part  of  the 
experimental  portion  of  this  contractual  effort.  Details  of  the 
experimental  work,  performed  by  the  Allison  Gas  Turbine  Division,  GMC,  are 
reported  in  Ref.  25.  In  addition  details  of  the  experimental  techniques, 
data  reduction,  etc.  can  be  found  in  the  above-mentioned  reference.  Both 
seal  configurations  considered  were  tested  at  pressure  ratios  of  2.0,  and 
these  are  the  cases  that  will  be  considered  in  this  report.  Schematics  of 
the  two  seals  are  shown  in  Figs.  16  and  17.  Probe  locations  are  noted  by 
the  stations  A,  B,  etc. 

For  the  straight-through  seal  configuration  probes  were  made  in  the 
centers  of  the  first  and  third  knife  tip  clearance  gaps,  i.e.,  stations  B 
and  I.  Probes  were  also  made  at  0.20  Inches  upstream  and  downstream  of  the 
edges  of  the  knife  tips,  i.e.,  at  stations  A,  C,  E,  G,  H,  and  J  • 
Additionally,  a  probe  was  made  at  the  half  point  of  the  first  cavity,  i.e., 
station  D  .  Measurements  in  the  clearance  gaps  consist  of  the  streamwise 
velocity  component  while  measurements  fore  and  aft  of  the  knife  tips  and  in 
the  first  cavity  consist  of  both  streamwise  and  transverse  velocity 
components. 

For  the  three  knife  stepped  seal  with  tapered  knives,  probes  were  made 
in  the  centers  of  each  knife  tip  In  the  clearance  gaps,  I.e.,  stations  B, 

F,  and  I  of  Fig.  17.  The  measured  results  for  this  seal  consist  of  the 
streamwise  velocity  component  only. 

The  calculations  used  for  the  comparison  were  previously  described  In 
the  section  on  the  calculation  of  leakage  rates,  specifically  cases  RC 
and  14B.  Both  of  the  computed  cases  utilized  a  k-e  turbulence  model  to 
account  for  the  effects  of  turbulence.  The  governing  equations  consisted 
of  the  two  momentum  equations,  the  continuity  equation  and  the  energy 
equation.  Boundary  conditions  were  as  previously  described.  Cases  RC 


and  14B  utilized  the  mixing  length  solutions  of  cases  8B  and  14A 
respectively  as  Initial  conditions.  The  two  turbulence  equations,  the 
turbulence  kinetic  energy  equation,  Eq»  (27),  and  the  dissipation  of 
turbulence  kinetic  energy  equations,  Eq.  (28),  were  then  solved  to 
convergence  with  the  fluid  dynamic  variables,  u,  w,  p  and  h  frozen  to 
obtain  initial  values  of  k  and  e.  Then  the  fluid  dynamic  and  the  k  and  e 
equations  were  simultaneously  solved  until  a  steady  state  solution  was 
obtained. 

The  results  for  the  three  knife  straight-through  seal  with  tapered 
knives  are  presented  in  Figs.  18  through  23.  For  the  flow  across  the 
centerline  of  the  clearance  gap  at  the  tip  of  the  first  knife  (Fig.  18), 
the  calculated  flow  is  uniformly  of  lower  velocity.  Qualitatively  the 
shape  of  the  streamwise  velocity  profile  is  similar  with  both  measured  and 
calculated  profiles  showing  a  tendency  of  the.  flow  to  separate  on  the  tip 
of  the  knife.  This  can  be  seen  by  the  shape  of  the  streamwise  velocity 
profile  on  the  knife  tip  having  a  significantly  smaller  gradient  (skin 
friction)  than  on  the  land.  The  measured  velocity  profile  shows  a  thinner 
boundary  layer  on  both  the  knife  tip  and  the  land  than  the  predicted 
boundary  layer  profile.  The  measured  profile  shows  no  discernable  boundary 
layer  on  the  land.  The  difference  in  the  magnitude  of  the  calculated  and 
measured  streamwise  velocity  is  consistent  with  the  results  of  the  leakage 
study  where  the  calculated  leakage  rate  was  lower  than  the  measured  value 
(as  was  previously  discussed).  The  results  for  the  streamwise  velocity 
profile  in  the  clearance  gap  at  the  third  knife  tip  (Fig.  19)  again  shows 
uniformly  lower  calculated  values  than  those  measured.  In  this  case  the 
measured  profile  shows  a  monotonlcaliy  Increasing  streamwise  velocity  as 
the  distance  from  the  knife  tip  is  Increased.  The  calculated  velocity 
profile  shows  a  slightly  skewed  profile  around  the  centerline  in  the  middle 
of  the  gap  with  the  tip  boundary  layer  being  somewhat  thicker  than  the  land 
boundary  layer.  Again  the  meaoured  results  show  no  discernable  boundary 
layer  on  the  land. 

The  results  in  the  region  0.20  Inches  upstream  of  the  leading  edge  of 
the  first  knife  are  presented  in  Fig.  20  and  the  tabulated  data  are 
presented  in  Tsble  3.  Qualitatively  the  experimental  data  agree  with  the 
calculated  flow  angles.  (If  the  calibration  of  the  hot-wire  Is  valid,  flow 
angles  and  velocity  magnitudes  should  be  accurate  to  a  few  percent). 
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The  velocity  magnitudes  calculated  are  uniformly  lower  than  those  measured 
which  is  again  consistent  with  the  calculated  leakage  rate  being  lower  than 
the  measured  values. 

The  results  in  the  cavity  regions  between  the  knives  are  presented  in 
Figs  21  and  22.  Fig.  21  shows  the  calculated  and  measured  velocity  vectors 
in  the  first  cavity  and  Fig  22  shows  the  calculated  and  measured  velocity 
vectors  in  the  second  cavity.  The  corresponding  tabular  values  are 
presented  in  Table  3.  In  the  cavity  between  the  first  and  second  knives 
the  calculated  and  measured  flow  angles  are  quite  similar,  viz.,  the  flow 
angles  very  nearly  match.  This  indicates  that  the  sizes  and  locations  of 
the  recirculation  zone  are  similar.  Generally  the  measured  values  of  the 
velocity  in  the  cavity  are  higher  than  the  calculated  values.  However,  in 
the  region  above  the  recirculation  zone  the  calculated  magnitude  of  the 
velocity  is  higher  at  stations  C  and  D  and  lower  at  station  E. 

In  the  cavity  between  the  second  and  third  knives  the  results  are 
similar  to  those  in  the  first  cavity.  The  flow  angles  predicted  by  the 
calculation  are  a  good  approximation  of  the  measured  values,  sec 
Table  3.  Again  the  measured  speeds  in  the  recirculation  zone  are  higher 
than  those  calculated.  The  behavior  of  the  speeds  above  the  recirculation 
zone  is  the  same  as  in  the  flrat  cavity,  l.e,,  at  station  G  the  calculated 
speeds  are  higher  than  measured  end  at  station  H  the  calculated  speeds  are 
lower  than  the  measured  speeds. 

At  the  last  station  where  speeds  and  flow  angles  were  measured, 
station  J  ,  measurements  show  a  considerably  different  flow  structure  than 
that  calculated  (see  Pig.  23  and  Table  3).  The  calculated  flow  angles  do 
not  In  general  agree  with  the  measured  values.  The  calculated  speeds 
behind  the  knife  are  of  approximately  Che  same  magnitude  as  those 
measured.  However  the  calculated  speed  at  the  upper  moat  data  point  is 
over  three  times  the  measured  value.  The  calculation  predicts  a  larger 
recirculation  zone  downstream  of  last  knife  than  the  measured  data  seem  to 
predict.  Although  It  is  difficult  to  determine  from  the  limited  amount  of 
data  available  in  this  region,  the  measured  recirculation  zone  appears  to 
be  much  thinner  than  the  predicted  one. 

The  results  for  STLD  flow  through  the  three  knife  stepped  seat  with 
tapered  knives  are  presented  In  Figs.  24  through  26.  Qualitatively  the 
computed  and  measured  values  of  the  screamwise  velocity  profiles  in  the 


clearance  gaps  above  the  first  and  second  knife  tips  are  similar.  In  the 
first  knife  tip  clearance  gap  (Pig.  24)  the  peak  velocities  are  of 
approximately  the  same  value,  while  in  the  second  knife  tip  clearance  gap 
the  calculated  peak  velocity  is  16  percent  lower  than  the  measured  value. 

At  the  first  knife  tip  no  discernible  land  boundary  layer  is  observed 
in  the  measurements.  In  the  case  of  the  stepped  seal  (which  is  contrary  to 
the  case  of  the  straight-through  seal)  the  calculated  flow  has  a  somewhat 
lesser  tendency  to  separate  over  the  leading  edge  of  the  first  knife  than 
does  the  measured  flow. 

For  the  second  knife  (Fig.  25)  the  calculated  boundary  layer  profile 
over  the  knife  tip  is  very  similar  to  the  measured  one. 

For  the  third  knife  tip  clearance  gap  (Fig.  26)  the  measured  flow  is 
qualitatively  similar  to  that  measured  for  the  previously  described 
straight-through  seal.  A  monotonlcally  increasing  streamwise  velocity 
profile  with  increasing  normal  distance  from  the  knife  tip  is  observed  in 
both  cases.  Again  no  discernible  measured  boundary  layer  is  seen  on  the 
land.  Although  qualitatively  the  calculated  and  measured  profiles  are 
dissimilar,  the  peak  speeds  for  the  two  profiles  are  within  a  few  percent 
of  each  other.  As  with  the  three  knife  6traight-through  seal  with  tapered 
knives,  the  measured  data  would  seem  to  indicate  a  smaller  recirculation 
tone  downstream  of  the  laat  stepped  seal  knife  and  hence  a  more  rapid 
expansion  than  that  predicted  by  the  calculation. 

In  general  the  comparisons  between  the  measured  and  the  calculated 
flows  are  encouraging.  Since  this  is  the  first  effort  to  both  predict  the 
flow  in  such  environments  and  to  perform  detailed  measurements  on  these 
types  of  labyrinth  seals,  the  results  are  better  than  might  Have  been 
anticipated.  More  experience  in  both  the  experimental  and  analytical 
efforts  should  lead  to  considerable  Improvements  in  both  areas.  An 
Important  point  to  be  mentioned  at  this  Juncture  Is  that  no  attempt  was 
made  to  'fine  tune'  the  calculations  to  get  better  comparisons  with  the 
experimental  data,  the  analysis  models  were  in  fact  run  before  the  test 
data  were  available,  the  measured  data  comparison  was  performed  with  model 
calculations  that  were  not  influenced  by  the  experimental  program. 

The  type  of  'tuning'  which  could  be  examined  would  concern  improving 
the  numerical  accuracy  and  modifying  the  Inflow  profiles  to  more  accurately 
reflect  the  conditions  (apparently)  occurring  in  the  experimental  seals. 


3.4  Rotating  Labyrinth  Seal  Calculation 


The  first  three  categories  of  calculations  performed  in  this  study 
were  concerned  with  "two-dimensional"  labyrinth  seal  configurations  having 
a  rectangular  clearance  gap.  All  experimental  cases  were  run  with  a 
spanwise  dimension  of  6.28  inches  to  simulate  the  static  test  rig. 

The  flow  through  the  clearance  gaps  would  essentially  be  two-dimensional  at 
the  experimental  aspect  ratios,  i.e.,  the  endwall  effects  were  neglected. 
Most  practical  applications  of  labyrinth  seals  are  for  rotating  equipment 
in  which  case  an  axially  symmetric  set  of  coordinates  must  be  used  to 
describe  the  seal  geometry.  In  order  to  demonstrate  this  capability  a 
sample  case  was  run  with  rotation. 

The  three  knife  straight  seal  with  tapered  knives  was  chosen  with  the 
knife  side  rotating  for  demonstration  at  a  pressure  ratio  of  2.0.  As  in 
the  nonrotating  case  the  flow  was  initially  assumed  to  be  stagnant  and  the 
back  pressure  lowered.  The  flow  was  then  drawn  through  the  seal  and  the 
basic  nonrotating  flow  pattern  established.  Then  the  rotor  was  turned 
until  the  desired  rotational  speed  was  obtained  and  the  converged  solution 
with  rotation  was  obtained.  The  equations  solved  are  the  transformed 
cylindrical  polar  Navier-Stokes  equations.  The  transformation  is  required 
because  of  the  tapered  knives  which  yield  a  nonorthogonal  coordinate 
system.  Three  separate  momentum  aquations,  the  continuity  equation  and 
an  energy  equation  are  solved.  Because  of  the  rotation  of  the  rotor, 
a  rotational  momentum  equation  must  be  solved  even  though  all  rotational 
derivatives  are  zero.  The  physical  dimensions  for  the  seal  configuration 
chosen  are  shown  in  Fig.  3g.  In  addition  the  radius  of  the  inner 
cylindrical  surface  of  the  rotor  was  arbitrarily  chosen  as  0.254  meters 
(or  100  clearances)  and  the  rotational  speed  was  chosen  as  6,000 
revolutions  per  minute,  628  radians  per  sec. 

A  sample  t,  -eamline  plot  in  the  vicinity  of  the  first  knife  for  the 
rotating  seal  case  is  shown  in  Fig.  27.  Qualitatively  the  streamline 
patterns  are  the  same  as  for  the  nonrotating  case.  The  steady 
state  value  of  the  flow  parameter  for  the  rotating  case  was 
0.315  Ibm  °rI /2/ibf-sec  V8  0.331  lbm  °R^/^ibf-sec  for  the  nonrotating 
case.  Thus  the  calculated  overall  effect  of  the  cylindrical  symmetry  and 
the  rotational  effects  is  to  decrease  the  value  of  the  flow  parameter 
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(and  hence  comparable  leakage  rate)  for  this  seal  configuration  by  about 
5  percent  over  the  nonrotating  case. 

Another  effect  predicted  by  the  calculation  procedure  was  the  rather 
large  amount  of  swirl  velocity  that  exists  far  downstream  of  the  last 
knife.  Since  the  height  of  the  knives  are  11/12  of  the  total  distance  from 
the  rotor  cylinder  to  the  land,  the  knives  impart  large  amounts  of  swirl 
into  the  upper  portion  of  the  flow  domain  (the  knife  tip  speeds  are  12.6 
times  the  mean  inlet  velocity  and  71  percent  of  the  peak  streamwise 
velocity  downstream  of  the  last  knife  tip).  In  the  immediate  vicinity  fo 
the  trailing  edge  of  the  last  knife,  the  swirl  is  rather  rapidly  dissipated 
by  the  large  viscous  effects  (at  one  clearance  downstream  the  swirl 
velocity  has  decreased  to  45  percent  of  the  peak  tip  velocity).  However 
downstream  of  this  region  the  losses  are  much  smaller.  At  100  clearances 
downstream  the  swirl  velocity  drops  by  18  percent.  Thereafter  the 
decreases  in  the  swirl  level  become  essentially  zero.  This  is  probably  due 
to  the  small  viscous  forces  in  the  core  of  the  flow  far  downstream  and  the 
lack  of  a  circumferential  pressure  gradient.  On  the  other  hand  in  the 
region  upstream  of  the  first  knife,  the  extent  of  the  penetration  of  the 
swirl  into  the  flow  is  minimal*  The  swirl  imparted  into  the  flow  by  the 
rotor  is  small  upstream  of  the  first  knife.  As  the  first  knife  is 
approached  some  of  the  swirl  is  convected  away  from  the  rotor  surface  by 
the  increasing  transverse  velocity  (caused  by  the  turning  of  the  streamwise 
velocity  to  flow  over  the  knife). 

The  general  observation  for  the  rotation  case  is  that  the  code  appears 
to  be  able  to  successfully  calculate  seal  flows  in  the  presence  of 
rotation.  In  future  work  this  should  probably  be  an  area  where 
considerable  effort  should  be  expended  both  experimentally  and 
computationally  as  this  is  the  actual  environment  in  which  real  seals 
exist. 


4.0  CONCLUSIONS 


The  results  of  this  first  effort  to  calculate  the  flow  In  labyrinth 
seals  by  the  numerical  solution  of  the  Navier-Stokes  equations  Is  very 
encouraging.  It  has  been  demonstrated  that  the  flow  In  a  wide  variety  of 
realistic  labyrinth  seal  geometries  can  be  calculated  under  various  flow 
conditions.  The  flow  In  both  straight-through  and  stepped  seal  geometries 
has  been  successfully  calculated.  Both  orthogonal  and  nonorthogonal 
coordinate  systems  have  been  used,  and  the  flow  has  been  calculated  for 
both  a  planar  and  axlsymmetric  system.  Pressure  ratios  as  large  as  five  to 
one  have  been  calculated  with  no  apparent  problem.  Seals  with  multiple 
knives  have  been  considered,  and  the  flow  was  successfully  predicted.  A 
variety  of  boundary  conditions  have  been  successfully  utilized  In  these 
calculations,  and  a  general  starting  procedure  has  been  developed  that  can 
be  used  with  any  seal  geometry  and  for  any  pressure  ratio.  The  numerical 
procedure  has  proven  to  be  robust,  l.e.,  all  calculations  that  were 
attempted  produced  converged  solutions.  Both  mixing  length  and 
two-equation  turbulence  models  were  successfully  used  for  the 
calculations.  Numerical  difficulties  often  associated  with  the 
two-equation  (k-e)  turbulence  model  were  eliminated  and  converged  solutions 
obtained  for  cases  that  had  not  previously  converged  with  the  use  of  the 
k-e  turbulence  model.  Calculation  of  the  flow  for  an  axlsymmetric  rotating 
labyrinth  seal  presented  no  problem  cither  with  the  physics  or  in 
numerically  converging  the  case. 

Comparison  of  the  calculated  results  with  experimental  results  was  In 
general  very  encouraging,  especially  when  It  is  remembered  that  this  is  the 
first  effort  for  these  classes  of  geometries.  The  computer  code  has 
demonstrated  an  ability  to  accurately  calculate  the  leakage  rates  for  a 
w'.de  variety  of  geometries  and  flow  conditions.  The  two  notable  exceptions 
are  the  three  knife  straight-through  seal  with  tapered  knives  and  three 
knife  stepped  seal  with  recangular  knives  at  pleasure  rat ’ os  of  2.0. 

The  comparison  of  the  hot-wire  data  with  the  calculated  results  for  the 
three  knife  straight-through  seal  with  tapered  knives  shows  qualitatively 
similar  results.  The  three  knife  stepped  seal  with  tapered  knives 
atmwn  good  >|tml  Itutlvi'  agreement  between  experiment  and  calculation. 

The  performance  curves  predicted  by  the  calculation  procedure  generated 
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curves  similar  to  experimental  data  and  to  those  predicted  by  the  Allison 
design  model.  A  major  advantage  of  the  calculation  procedure  would  be  its 
ability  to  calculate  these  performance  curves  for  labyrinth  seals  for  which 
there  is  no  data  base. 

The  calculations  performed  under  this  effort  required  a  reasonable 
amount  of  computer  time*  Further,  it  is  expected  that  in  the  near  future 
with  a  vectorized  program  typical  calculations  could  be  run  for  on  the 
order  of  less  than  $200  per  data  point.  It  is  hoped  that  the  computer  code 
will  be  integrated  into  the  design  process  in  the  near  future.  Because  of 
the  very  general  nature  of  the  computer  code,  i.e.,  it  has  the  ability  to 
accept  any  reasonable  coordinate  system  and  to  perform  calculations  in  that 
system,  it  would  be  desirable  to  utilize  the  code  to  perform  calculations 
for  advanced  seal  concepts.  Finally,  it  would  be  desirable  to  perform 
calculations  for  actual  full-scale  labyrinth  seals  in  the  rotating  mode  as 
this  is  the  environment  that  is  actually  experienced  in  sealing 
applications. 
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LIST  OF  SYMBOLS 


Symbols 

A 

CP 

C1 

c2 

Cy 

D 

F 

H 
h 
J 
k 
L 

*m 

a 

M 

+ 
n 

+ 

n8 
P 

Pr 

R 
R 
Re 
r 
S 


Area  or  matrix  of  time  linearization  coefficients 

Specific  heat  at  constant  pressure 

Constant  for  two  equation  turbulence  model 

Constant  for  two  equation  turbulence  model 

Constant  for  two  equation  turbulence  model 

Rate  of  deformation  tensor  or  elements  of  spatial 
differential  operators 

Vector  of  convection  and  diffusion  terms 
<Eqs.<3l)-<33)) 

Vector  of  time  terms 

Enthalpy 

Jacobian 

Turbulence  kinetic  energy 

Matrix  of  linear  differential  operators 

Mixing  length 

Mass  flux 

Mixture  molecular  weight 

Unit  vector  in  normal  direction 

Unit  vector  in  symmetry  direction 

Pressure 

Prandtl  number 

Heat  flux  vector 

Universal  gas  constant 

Reynolds  number 

Seal  expansion  ratio 

Vector  of  source  terms  (Eq*(34)) 
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XI  s  <+ 


Symbols 


t  Time 

T  Temperature 

U  Velocity  component 

A 

U  Specific  internal  energy 

u  Cartesian  velocity  component 

V  Average  velocity 

Velocity  vector 

Vector  of  flux  variables  (Eq.(30))  or  mass  flux 

Cartesian  coordinate 
y  General  coordinate 


Greak  Symbols 
6 
A 
fi 
c 
0 
•c 

U 

V 

0 

0 

t 

♦ 

4 


Crank-Nicolson  factor 
Change 

Kronecker  delta 

Dissipation  of  turbulence  kinetic  energy 

Rotational  direction 

Thermal  conductivity 

Dynamic  viscosity 

Kinematic  viscosity 

Density 

Dissipation  parameter 
Stress  censor 
Rnerity  dissipation 

Vector  of  dependent  variables  or  flow  parameter 
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Subscripts 


art 

call 

corr 

D 

£ 

n 

o 

s 

T 

Tl 

T2 

w 


Artificial 

Calculated 

Correlated 

Downstream 

Laminar 

Normal 

Stagnation 

Symmetry 

Turbulent 

First  tangential  direction 
Second  tangential  direction 
Wall 


*i 

*3 

1 

2 
3 


Associated  with  ifch 

Associated  with 
First  direction 
Second  direction 
Third  direction 


Cartesian  direction 
Cartesian  direction 


Superscripts 

n 

T 

* 

** 


1 

2 

3 


nth  time  step 
Transpose 

First  intermediate  time  level 

Second  intermediate  time  level 

Fluctuation 

First  direction 
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STEPPED  SEALS 


STRAIGHT  SEALS 

CL 


FiRurc  2  -  Labyrinth  Seal  Nomenclature. 
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Figure  3  -  Straight -Through  Seal  Configurations. 


0.1 


i 


0.1 


(ft)  Tapered  -  Throe  Kttlvey 
Figure  3  -  (Continued) 
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Figure  6  -  Pressure  Comparisons  -  Three  Knife  Straight-Through  Seal  with  Tapered  Knives. 


A  DATA 

•  CALCULATION 


Velocity  Vectors  -  Single  Knife  Straight-Through  Seal  with 
Tapered  Knife. 


Figure  9  -  Streamlines  -  Three  Knife  Straight-Through  Seal  with 
Tapered  Knives. 


Figure  10  -  Mach  Number  Contours  -  Single  Knife  Straight-Through  Seal 
with  a  Slanted  Knife. 
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Figure  14  -  Performance  Plot  -  Three  knife  Straight-Through  Seal  With  Tapered  Knives. 


.25 


Figure  16  -  Data  Stations  -  Three  Knife  Straight-Through  Seal 
with  Tapered  Knives. 


ca  Stations  -  Three  Knife  Stepped  Seal  with  Tapered  Knives 


Figure  19  -  Streamvise  Velocity  Comparison  -  Station  I  -  Three  Knife 
Straight-Through  Seal  with  Tapered  Knives. 
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VELOCITY  SCALE:  l"=  20  m/scc 

-o  20-  Velocity  Component  Comparison  -  Station  A  -  Three  Knife 
St  ra  H>,ht  -Through  Seal  With  Tapered  Knives  . 
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VELOCITY  SCALE  t  i  40  m/wc 

,:iguro  21  -  Velocity  Component  Comparison  -  First  Cavity  (St  at  ions  C ,  I) ,  anil  I'.) 
Throe  Knife  St  ra  iglit  -Tlirout'h  Seal  with  Tapered  Knives. 
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POSITION  -  INCHES  UP  FROM  FLOOR 


0.4  < — 


VELOCITY  SCALE:  I  =  40  m/sec 

22  "  comparison  -  So, -on, I  Cavity  (Stations  G  and  H) 

lireo  Knife  StraiRht-lhrouRh  Seal  With  Tapered  Knives. 


Figure  25  -  Streamwise  Velocity  Comparison  -  Station  F  -  Three  Knife 
Stepped  Seal  with  Tapered  Knives. 


Scralgbc  Through  Labyrinth  Seals 


ocity  Component  Comparison  -  Three  Knife  Straight-Through 
1  with  Tapered  Knives  sc  P  /P_«2.G. 
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APPENDIX  A  -  TRANSFORMATION 


It  is  desired  to  show  that  given  the  general  transformation 

y*  *  y* <*,,*,» *s)  (A-i) 

that  the  relationship 

-^r  (d  Ill)  •  0  (A-2) 

dyJ  d*{ 

is  valid*  Where  m  in  Section  2.2,  J  is  the  Jacobian  of  the  transformation 
Kq»  (A-I).  Applying  the  chain  rule  to  8q,  (A-i)  yields  the  relationships 


dy1  4y*  41, 

4,'  *  11  '  «„  4y~ 


(A-3) 


Writing  Eq.  (A-3)  for  each  of  the  three  k  directions  for  one  yJ  yields 
three  linear  relationships  for  the  unknowns  dyi/Sx^.  These  can  bo 
solved  by  applying  Cramers  rule  to  yield  for  j-1,2  and  3  respectively 
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d/ 


ay1 


«* 

dy* 

dy3  dy2 

dfs 

dy* 

dy3  * 

dy3 
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ay2 

a*2  a?3 

ax2  a*5 

(A-7) 

dx, 

ay5  ay* 

dy*  dy5 

ay2 

—  g 

ax3  ax, 

dx5  a?, 

a*2 

ay5  dy‘ 

dy*  dy3 

(A-8) 

ay2 

dx,  dx2 

dx,  dx2 

a*3 

dy3  ay1 

ay*  dy3 

(A-9) 

ay*  ^ 

a* g  axa 

dx2  ax3 

(A-10) 
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ay*  dy*  " 
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a*3  ai, 

ax3  a?, 

(A-li) 

a*2 

ay*  dy2 

dy2  dy* 

ay3 

ax,  dxj. 

dx,  dx2 

(A- 12) 

dy*  dy2 

dy2  dy* 

Substitution  of  relationships  (A-4)  through  C A— 12)  into  Eq.  (A-2) 
substantiates  the  validity  of  Eq.  (A-2). 
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These  relationships  (A-4)  through  (A-12)  and  the  definition  of  the 
Jacobian,  J,  (Eq.  42),  are  used  to  calculate  the  geometric  groupings  that 
occur  in  the  governing  system  of  coupled  partial  differential  equations 
represented  by  Eq.  (46)  and  the  auxiliary  relationships  represented  by 
Eqs.  (47)~(50).  These  geometric  groupings  are  calculated  by  finite 
difference  techniques  in  an  analogous  manner  to  that  used  for  the  fluid 
dynamic  derivatives.  A  finite  difference  grid  distribution  is  setup  in  the 
computational  domain  (yj)  and  grid  points  are  associated  with  the 
cartesian  location  (xj).  All  derivatives  of  the  form  3xi/3yj  are 
then  approximated  by  their  finite  difference  analog,  i.e.,  central 
differences  for  interior  grid  points  and  three  point  one-sided 
approximations  on  boundaries  Once  these  derivatives  are  calculated  for 
all  values  of  i  and  j  the  Jacobian  J  and  the  geometric  groupings, 
represented  by  Eq.  (A-4)  through  Eq.  (A-12)  can  be  calculated. 
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APPENDIX  B  -  SOLUTION  PROCEDURE 
Background 

The  solution  procedure  employs  a  consistently-split  linearized  block 
implicit  (LBI)  algorithm  which  has  been  discussed  in  detail  in  Refs.  12 
and  26.  There  are  two  important  elements  of  this  method: 

(1)  the  use  of  a  noniterative  formal  time  linearization  to 
produce  a  fully-coupled  linear  multidimensional  scheme  which 
is  written  in  "block  implicit"  form;  and 

(2)  solution  of  this  linearized  coupled  scheme  using  a  consistent 
"splitting"  (ADI  scheme)  patterned  after  the  Douglas-Gunn  (Ref.  27) 
treatment  of  scalar  ADI  schemes. 

The  method  is  thus  referred  to  as  a  split  linearized  block  implicit  (LBI) 
scheme.  The  method  has  several  attributes: 

(1)  the  noniterative  linearization  is  efficient; 

(2)  the  fully-coupled  linearized  algorithm  eliminates  instabilities 
and/or  extremely  slow  convergence  rates  often  attributed  to  methods 
which  employ  ad  hoc  decoupling  and  linearization  assumptions  to 
identify  nonlinear  coefficients  which  are  then  treated  by  lag  and 
update  techniques; 

(3)  the  splitting  or  ADI  technique  produces  an  efficient  algorithm 
which  is  stable  for  large  time  steps  and  also  provides  a  means  for 
convergence  acceleration  for  further  efficiency  in  computing  steady 
solutions; 

(4)  intermediate  steps  of  the  splitting  are  consistent  with  the 
governing  equations,  and  this  means  that  the  "physical"  boundary 
conditions  can  be  used  for  the  intermediate  solutions.  Other 
splittings  which  are  inconsistent  can  have  several  difficulties  in 
satisfying  physical  boundary  conditions  (Ref.  12). 

(3)  the  convergence  rate  and  overall  efficiency  of  the  algorithm  are 
much  less  sensitive  to  mesh  refinement  and  redistribution  than 
algorithms  based  on  explicit  schemes  or  which  employ  ad  hoc 
decoupling  and  linearization  assumptions.  This  Is  Important  for 
accuracy  and  for  computing  turbulent  flows  with  viscous  sublayer 
resolution;  and 
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(6)  the  method  is  general  and  is  specifically  designed  for  the 

complex  systems  of  equations  which  govern  multiscale  viscous  flow 
in  complicated  geometries. 

This  same  algorithm  was  later  considered  by  Beam  and  Warming  (Ref.  25),  but 
the  ADI  splitting  was  derived  by  approximate  factorization  instead  of  the 
Douglas-Gunn  procedure.  They  refer  to  the  algorithm  as  a  "delta  form" 
approximate  factorization  scheme.  This  scheme  replaced  an  earlier  non-delta 
form  scheme  (Ref.  34)  which  has  inconsistent  intermediate  steps. 

Split  LBI  Algorithm 
Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  three  or  four 
equations:  continuity  and  two  or  three  components  of  the  momentum  equation 
in  three  or  four  dependent  variables:  p,  u,  v,  and/or  w.  Using  notation 
similar  to  that  in  (Ref.  12),  at  a  single  grid  point  this  system  of  equations 
can  be  written  in  the  following  form: 

3H(*)/3t  -  D($)  +  S(*)  (B-3) 

where  $  is  the  column-vector  of  dependent  variables,  H  and  S  are  column- 
vector  algebraic  functions  of  ♦,  and  D  is  a  column  vector  whose  elements  are 
the  spatial  differential  operators  which  generate  all  spatial  derivatives 
appearing  in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  basod  on  the  following  two-level  Implicit 
time-difference  approximations  of  (B-3): 

Hn)/4t  -  B(Dn*l+  Sn+l)  ♦  ( 1-8)  (Dn  +  S")  (8-4) 

where,  for  example,  Hn+*  denotes  M(4n'*'l)  and  At  *  tn+*  -  tn.  The 
parameter  B  (0.5  <  8  £  1)  permits  a  variable  time-centering  of  the  scheme, 
with  a  truncation  error  of  order  jAt**,  (8  -  1/2)  At). 

A  local  time  linearization  (Taylor  expansion  about  ♦n)  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear  dlfferon- 
ttal  operator  L  (cf.  Ref.  12)  such  that 

Dn*1  -  n"  +  inun+1-  ♦”)  ♦  0(At2)  (B-5) 
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Similarly 


Hn+I  «  Hn+  (3H/3^)n  <$n+1  -*")  +  <)  (At2)  (B-6) 

Sn+1  -  Sn+  (3S/3$)n  ($n+1  -  $n)  +  0  (At2)  (B-7) 

Eqs.  (B-5  through  B-7)  are  inserted  into  Eq.  (B-4)  to  obtain  the  following 
system  which  is  linear  in 

(A  -  BAt  L°)  ($°+1  -  $n)  »  At  (Dn  +  S")  (B-8) 

and  which  is  termed  a  linearized  block  implilcit  (LBI)  scheme.  Here,  A 
denotes  a  matrix  defined  by 

A  =  (3H/3$)n  -  BAt  (3S/3<»n  (B-9) 


Eq.  (B-8)  has  0  (At)  accuracy  unless  H  =  in  which  case  the  accuracy  is  the 
same  as  Eq.  (B-4). 

Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate 
cross-derivative  terms  and  also  turbulent  viscosity  and  artificial  dissipa¬ 
tion  coefficients  which  depend  on  the  solution  variables.  Although  formal 
linearization  of  the  convection  and  pressure  gradient  terms  and  the  resulting 
implicit  coupling  of  variables  is  critical  to  the  stability  and  rapid  con¬ 
vergence  of  the  algorithm,  this  does  not  appear  to  be  Important  for  the 
turbulent  viscosity  and  artificial  dissipation  coefficients.  Since  the 
relationship  between  ue  and  dj  and  the  mean  flow  variables  is  not  conven¬ 
iently  linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at 
tn  during  each  time  step,  Notat tonally,  this  is  equivalent  to  neglecting 
terms  proportional  to  3  uc/3#  or  3dj/3$  In  Un,  which  are  formally  pre¬ 
sent  In  the  Taylor  expansion  (8-5),  while  retaining  all  terms  proportional  to 
u0  or  dj  In  both  Ln  and  Dn. 

It  has  been  found  through  extensive  experience  that  this  has  little  if 
any  effect  on  the  performance  of  the  algorithm.  This  treatment  also  has  the 
added  benefit  that  the  turbulence  model  equations  can  be  decoupled  from  the 
system  of  mean  flow  equations  by  an  appropriate  matrix  partitioning 
(cf.  Ref.  26)  and  solved  separately  ln  each  step  of  the  ADI  solution 
procedure.  This  reduces  the  block  size  of  the  block  trldlagonal  systems 
which  must  be  solved  ln  each  step  and  thus  reduces  the  computational  labor. 

In  addition,  the  viscous  terms  ln  the  present  formulation  Include 
a  number  of  cross-derivative  terms  Implicitly  within  the  ADI  treatment 
which  follows.  It  is  not  at  all  convenient  to  handle  these  implicit 
cross-derivative  terms;  and  consequently,  all  cross-derivative 


terms  are  evaluated  explicitly  at  tn.  For  a  scalar  model  equation  representing 
combined  convection  and  diffusion,  it  has  been  shown  by  Beam  and  Warming  that  the 
explicit  treatment  of  cross-derivative  terras  does  not  degrade  the  unconditional 
stability  of  the  present  algorithm.  To  preserve  notational  simplicity,  it  is 
understood  that  all  cross-derivative  terras  appearing  in  Ln  are  neglected  but 
are  retained  in  Dn.  It  is  important  to  note  that  neglecting  terms  in  Ln  has 
no  effect  on  steady  solutions  of  Eq.  (B-8),  since  =  0  and  thus  Eq.  (B-8) 

reduces  to  the  steady  form  of  the  equations:  Dn  +  Sn  =  0.  Aside  from 
stability  considerations,  the  only  effect  of  neglecting  terras  in  Ln  is  to 
introduce  an  0  (At)  truncation  error. 


Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (B-8)  is  split  using 
ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional  operator  L  is 
rewritten  as  the  sum  of  three  “one-dimensional**  sub-operators  Lj  (i  =  1,  2,  3) 
each  of  which  contains  all  terms  having  derivatives  with  respect  to  the  i-th 
coordinate.  The  split  form  of  Eq.  (B-8)  can  be  derived  either  by  following  the 
procedure  described  by  Douglas  and  Gunn  (Ref.  27)  in  their  generalization  and 
unification  of  scalar  ADI  schemes  (as  done  in  Refs.  12  and  26),  or  by  using 
approximate  factorization.  For  the  present  system  of  equations,  the  split 
algorithm  is  given  by 


(A  -  SAthj)  (* 
(A  -  8AtLf2')  (♦ 
(A  -  flAthj)  <$ 
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♦">  ■ 
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(B-lOa) 

** 

-  ♦"> 
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> 

♦ 
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IK 
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■  A  ($ 
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where  and  are  consistent  intermediate  solutions.  If  spatial 
derivatives  appearing  in  l.j  and  0  are  replace  by  three-point  difference 
formulas,  as  Indicated  previously,  then  each  step  ln  Eqs.  (B-10a,b  and  c)  can  he 
solved  by  a  block-tridiagonal  elimination. 

Combining  Kqs.  (B-10a,b  and  c)  gives 

(A  -  BAtl")  A~‘  (A  -  BAthj)  < A  -  BAtl.")  (♦’>+1  -  *”)  -  6t  <D°  +  S°)  (B-ll) 

which  approximates  the  unspllt  scheme,  Eq.  (B-8)  to  0  (At^).  Since  the 
intermediate  steps  are  also  consistent  approximations  for  Eq.  (B-8),  physical 
boundary  conditions  can  be  usod  for  and  (Refs.  12  and  26).  Finally, 
since  the  L|  are  homogeneous  operators,  it  follows  from  Eqs.  (B-I0a,b  and  c) 
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that  steady  solutions  have  the  property  that  «*  ■  $**  =*  and 

satisfy 


Dn  +  Sn  *  0 


<  B— 12) 


The  steady  solution  thus  depends  only  on  the  spatial  difference  approximations 
used  for  Eq.  (B-12),  and  does  not  depend  on  the  solution  algorithm  itself. 
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